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ABSTRACT 

Feedback from massive stars is believed to play a critical role in shaping the galaxy mass 
function, the structure of the interstellar medium (ISM), and the low efficiency of star for- 
mation, but the exact form of the feedback is uncertain. In this paper, the first in a series, we 
present and test a novel numerical implementation of stellar feedback resulting from momen- 
tum imparted to the ISM by radiation, supernovae, and stellar winds. We employ a realistic 
cooling function, and find that a large fraction of the gas cools to < 100 K, so that the ISM 
becomes highly inhomogeneous. Despite this, our simulated galaxies reach an approximate 
steady state, in which gas gravitationally collapses to form giant 'molecular' clouds (GMCs), 
dense clumps, and stars; subsequently, stellar feedback disperses the GMCs, repopulating 
the diffuse ISM. This collapse and dispersal cycle is seen in models of SMC-like dwarfs, 
the Milky-Way, and z ^ 2 clumpy disk analogues. The simulated global star formation ef- 
ficiencies are consistent with the observed Kennicutt-Schmidt relation. Moreover, the star 
formation rates are nearly independent of the numerically imposed high-density star forma- 
tion efficiency, density threshold, and density scaling. This is a consequence of the fact that, 
in our simulations, star formation is regulated by stellar feedback limiting the amount of very 
dense gas available for forming stars. In contrast, in simulations without stellar feedback, i.e., 
under the action of only gravity and gravitationally-induced turbulence, the ISM experiences 
runaway collapse to very high densities. In these simulations without feedback, the global star 
formation rates exceed observed galactic star formation rates by 1 — 2 orders of magnitude, 
demonstrating that stellar feedback is crucial to the regulation of star formation in galaxies. 

Key words: galaxies: formation — galaxies; evolution — star formation: general — cosmol- 
ogy; theory 



1 INTRODUCTION 

Feedback from massive stars plays a critical role in the evolution of 
galaxies. Cosmological models of galaxy evolution generically find 
that, without strong stellar feedback, the net stellar mass formed 
from cooled baryons exceeds that observed by an order of magni- 



tude or more, particularly in lower mass halos (e.g. Katz et al.|1996| 



Somerville & Primack'1999[ |Cole et al.|2000[|Springel & Hernquist| 
2003b ; Keres et al. 2009 and references therein). Related problems 
exist on smaller scales within galaxies. The observed relationship 
between star formation rate density and gas surface density - the 
Kennicutt-Schmidt (KS) law - implies that star formation is slow 
averaged over galaxies as a whole: the gas consumption timescale 
is ~ 50 dynamical times l |Kennicutt^l998) , much longer than the 
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naive estimate of ~ a few dynamical times one might expect in 
self-gravitating gas. Similar gas consumption times are found even 
in dense regions in galaxies (e.g., Krumholz & Tan 2007; but see 
also |Murray |20TT||Schruba et al. |20 1 0,.Feldmann & Gnedin]20TT}. 
Moreover, observations in the Milky Way and nearby galaxies show 
that individual giant molecular clouds (GMCs) - the sites of star 
formation - convert only a few percent of their mass into stars dur- 
ing their lifetimes ( Zuckerman & Evans 1974) [Williams & McKee| 
[T997{|Evans|1999||Evans et al.,2009b^ . One of the leading expla- 
nations for this low star formation efficiency is that stellar feedback 
disrupts GMCs once enough stars have formed. 

Numerical simulations of isolated galaxies and galaxy merg- 
ers, as well as cosmological "zoom-in" simulations of individual 
halos, can now reach the resolution required to resolve the for- 
mation of GMCs, ~ 1 - 100 pc (see e.g. |Boum aud et al.| |2010[ 

_ . „ _ ^^^^ 



[Saitoh et al.|[2008l [Dobbs et al.[|20TT] [Tasker & Tan,|2009) (note 
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that GMCs in massive gas-rich galaxies are ~ kpc in size, sig- 
nificantly larger than in the Milky Way). If simulations do not, 
however, include physics that disrupts GMCs, they do not have a 
physically self-consistent model of the interstellar medium (ISM) 
on such scales. All of the gas will be unrealistically locked up in 
dense gaseous/stellar clusters, instead of being recycled back into 
the more diffuse ISM. Given resolution limitations, most recipes in 
galaxy and cosmological-scale simulations have been developed to 
treat star formation and feedback in a restricted "sub-grid" manner. 
However, without more detailed models of this physics, it is dif- 
ficult to assess how appropriate the sub-grid prescriptions are for 
different galaxy types. Moreover, the assumptions of such models 
break down and are no longer meaningful at the spatial (<pc) or 
time evolution (< Myr) scales of individual GMCs and ISM sub- 
structure. In particular, whenever a numerical simulation has the 
resolution to resolve Ihs formation of bound gaseous structures like 
GMCs, we believe that it is equally critical to include physics that 
can potentially disrupt such GMCs. 

Protostellar jets, HII regions, stellar winds, radiation pressure 
from young stars, and supemovae all appear to be important sources 
of feedback and turbulence in the ISM of galaxies. In regions of low 
mass star formation it is likely that protostellar jets dominate, but 
for the ISM as a whole massive stars are the most important sources 
of feedback. In relatively low density gas, heating by photoioniza- 
tion, stellar winds, and supernovae is critical ( [McKee & Ostriker] 
1 1 977[ [Matzner|2002 For denser gas, however, which often corre- 
sponds to most of the mass in a galaxy, energy deposition is ineffec- 
tive; the cooling time {tcooI = kT/An ^ 3000(r/10"*A:)(lcm"V") 
yrs, where A ~ 10~^' ergcm^ s~' is the cooling function) is short 
compared to the dynamical time for all but the lowest density gas, 
so the energy deposited into the gas by stellar feedback is rapidly 
radiated away. Even in the Milky Way, the hot ISM contributes only 
~ 10% to the total ISM pressure jBoulares & Cox|1990y In con- 
trast, the momentum supplied by stellar luminosity, winds, and su- 
pernovae cannot be radiated away, and is the most important source 
of feedback for dense gas in galaxies (e.g., [Murray et al.|2010) . 

Although it is widely believed that stellar feedback is critical 
for understanding the self-regulation of star formation within galax- 
ies, and for the cosmological evolution of galaxies themselves, it 
is quite challenging to treat this in galaxy-scale simulations, espe- 
cially with the computational limitations faced by previous gener- 
ations of simulations. As a result, many simulations have made the 
problem tractable by adopting effective equation of state models 
in which feedback processes are treated implicitly (e.g. |Springel| 
|& Hemquist|2003a| [Teyssier et al.|2010} , accounting for the (un- 
resolved) multi-phase turbulent structure of the ISM with an "ef- 
fective" high sound speed. Unfortunately, in this case many of the 
net effects of stellar feedback are then put in by hand - one cannot 
predict, e.g., either how efficient feedback is in different systems 
or whether stellar feedback drives galactic winds. More broadly, 
without simulations that explicitly model feedback, it is difficult to 
evaluate the accuracy of the various subgrid treatments in the liter- 
ature. 

Galactic-scale simulations that do include stellar feedback 
explicitly have often been forced to alter the physics in signifi- 
cant ways in order to obtain a desired result. The most common 
treatment is to only include thermal gas heating from supemovae. 
However, thermal feedback is very inefficient in the dense regions 
where star formation occurs, and in the ISM more broadly in star- 
bursts and gas-rich high-redshift galaxies. These problems are com- 
pounded when simulations cannot resolve the ISM phase structure, 
and smooth together dense GMCs and diffuse gas into a single 



average density (greatly increasing/decreasing the cooling time in 
dense/diffuse gas, respectively). For this reason, in order to make 
supernova feedback have any significant effect (even in MW-like 
galaxies), simulators often "turn off" cooling (often along with star 
formation and/or other hydrodynamic processes) for an extended 
period of time, much longer than TcooI (cooling is typically sup- 
pressed for ~ lO' — 10** yr, i.e., for a time comparable to a galaxy 
dynamical time and ~ lO"* times longer than the actual cooling time 
at the same density; see e.g. IThacker & Couchman|2000] |Gover-| 



|nato et al.|2007[ [Brook et al.|20II| l. Other models explicitly dis 
able certain interactions between gas flagged as "cold" and "hot" 
or deposit feedback energy in a non-cooling reservoir that serves to 
move gas from cold to hot phases ( Scannapieco et alT[2008) . Even 
with these adjustments, many such models have found it difficult 
to drive winds and suppress star formation at the level needed to 
explain the galaxy mass function (especially at low masses) and 
observed star formation efficiencies (see e.g. Guo et al. 2010, Pow^ 
[ell et al.|20TT1 [Brook et al.|201 I|[Nagamine,2010| and references 
therein). 

Simulations with supernova feedback that do not "turn off" 
cooling have found that galactic outflows can only be driven if addi- 
tional physics is included. For example, [Ceverino & K lypin (2009| 
were able to drive galactic winds by requiring that supernovae ex- 
plode well outside of the GMCs in which they formed. However 
(as they acknowledge), although this may well be important for 
galactic winds, it leaves the problem of locally preventing runaway 
collapse of dense star forming regions. 

The inefficiency of supernova heating in dense gas is physi- 
cally correct.' It is thus by no means clear that turning off cooling 
is an appropriate resolution of the 'problem' that supernova rem- 
nants cool! Instead, we believe that this points to the importance of 
including the momentum supplied by stellar feedback processes. 
This momentum input into the ISM can drive strong turbulence 
and can itself contribute to unbinding gas from galaxies, even in 
the limit of very rapid cooling ( Murray et al.|2005] l. To date, how- 
ever, this has only been treated in a phenomenological way, given 
the limited resolution of previous simulations. In particular, in a 
widely used implementation in the Gadget SPH code, gas parti- 
cles are "kicked" into a "wind" at a rate proportional to either the 
star formation rate or young stellar mass, with the wind velocity 
set by hand as a constant or a multiple of the galaxy escape ve- 
locityjjSpnngel & Hemq uist|2003a[ Oppenhei mer & Dave|2008{ 
[Sales et al.[2010(|Genel et al.[2010[ [Dalla VecchiaTs 'chaye[2008| >. 
All hydrodynamic interactions (e.g., shocks and pressure forces) 
for the "wind particles" are turned off until they escape the galaxy 
(and when this is not done, the effects of winds are substantially 
suppressed; [Sales et ^[2010^ . This model is useful for studying 
the impact of galactic outflows on the intergalactic medium and 
the galaxy mass function but it is clearly limited (especially within 
galaxies) and cannot predict the nature and origin of these winds. 

Motivated by these considerations, this is the first in a series 
of papers studying stellar feedback in numerical models of galax- 
ies and the resulting implications for problems such as the origin of 
galactic winds, the physics of gas inflow in galaxy mergers, and the 
properties of the ISM in high redshift galaxies.^ Ultimately, we will 



' It is, of course, true that simulations do not resolve the full multi-phase 
ISM into which supernovae propagate and that this can enable supernova 

energy to propagate to larger distances. 

^ Movies of these simulations are available at [https : / / www . cf a ■ [ 

^harvard . edu/~phopkins /Site/Research . html| 



© 0000 RAS, MNRAS 000, 000-000 



Self-Regulated Star Formation in Galaxies 3 



present results that include simple models of supernova heating, 
HII regions, and radiation pressure from massive stars (produced 
by the absorption and scattering of UV and IR radiation on dust). 
Feedback from a central active galactic nucleus may also be impor- 
tant for understanding some aspects of star formation in galaxies 
- particularly the cessation of star formation in massive galaxies - 
but this is a separate problem that we do not consider in this paper. 

It is still not currently computationally feasible to include all 
of the physics of stellar feedback in simulations that focus on galac- 
tic scales. The methods we develop therefore still rely on sub-grid 
models, but at the sub-cluster or sub-GMC scale, as opposed to the 
galaxy scale. The fact that different feedback processes dominate 
under different physical conditions (e.g., density) highlights the im- 
portance of including a range of physical processes when studying 
the effects of stellar feedback on galaxies and galaxy formation. 
Nonetheless, in this paper, we focus on isolated (non-cosmological) 
galaxies and only include feedback by momentum deposition from 
massive stars. Our motivation for doing so is several-fold. First, 
our model for momentum-deposition is sufficiently different from 
existing treatments of stellar feedback in the literature that it re- 
quires a detailed explanation. More importantly, however, we show 
that this simple model is, by itself, able to explain the Kennicutt- 
Schmidt relation and the low star formation efficiency in galaxies. 
Moreover, the star formation rates in our model galaxies typically 
change by less than a factor of --^ 2 when we include additional 
feedback processes (though other properties of the galaxies can 
change substantially, such as the morphology and phase-structure 
of the ISM - this is particularly true for low mass galaxy models). 

The remainder of this paper is organized as follows. In ^we 
describe our method of implementing feedback due to the injec- 
tion of momentum by young, massive stars. The Appendix contains 
tests varying some of the parameters of our numerical method. The 
galaxy models we study are described in Table[T]and ^2.3| We then 
summarize our key results on the star formation histories and struc- 
tural properties of our model galaxies (^. In ^we show that these 
results do not depend strongly on the physics of star formation at 
high densities, the uncertain feedback parameters, and numerical 
resolution. We then show that our model galaxies are consistent 
with the observed Kennicutt-Schmidt relation (^. Finally, in ^ 
we summarize our results and discuss their implications. 

2 METHODOLOGY 

The methods we present are general and can be implemented in 
both Eulerian and Lagrangian simulations. The specific simula- 
tions we carried out were performed with the parallel TreeSPH 
code GADGET-3 l |Springel|2005| >, based on a conservative formula- 
tion of smoothed particle hydrodynamics (SPH), which conserves 
energy and entropy simultaneously even when smoothing lengths 
evolve adaptively (see e.g., |Springel & Hemqu ist 2002, Hernquist 
|I993[ [O'Shea et al.|2005). Th e detailed numerical methodology is 
described in [Springel j2003] >, |Springel & Hernquist| p003a[ l, and 
|Springel et ak (f2005'). Our simulations include stars, dark matter, 
and gas, with new implementations of stellar feedback; we describe 
the salient features of this additional physics below. These calcula- 
tions do not include models of black hole growth and feedback. 



by fine structure lines. Specifically we tabulate the cooling function 
A(r) from 1 - 10''K with CLOUDY, for a medium with density 
n = Icm^^, solar abundances, and with an ionizing background 
matching that at z = 0.^ This is similar to the approach in a number 
of other simulations (see e.g. 'Robertson & Kravtsov||2008[ |Wise| 
|& Ab el 2008 ; Ceverino & Klypin 2009 ) and gives identical results 
to the tabulated A(r) presented in |Sanchez-Salcedo et al. (2002i. 
We are not attempting to follow the ISM chemistry and thus ignore 
the dependence of the cooling on abundance and radiation field. For 
our problems of interest, the cooling rates even at these low temper- 
atures are uniformly much shorter than the dynamical times in all 
the systems we model; therefore, even large (factor ~ 5) changes in 
the cooling curve make no significant difference to our conclusions 
(we have checked this explicitly). 

Because we allow cooling to very low temperatures, we also 
must account for finite simulation resolution by including a pres- 
sure floor to prevent artificial numerical fragmentation when the 
Jeans mass is not resolved jTruelove et aL 1997| l. We adopt the pre- 



scription in [Robertson & Kravtsov (2008 i, which ensures that the 



Jeans length is always resolved with A'jcans smoothing lengths. This 
density-dependent pressure floor is 



Pjeans ~ 1 . 2 N^'^^^^^ J ' G /iLl 



(1) 



where 7 = 5/3, p is the local density, and hsm\ the smoothing 
length. We typically adopt Means ~ 10, but have experimented 
with Means =4—15 and find similar results (the morphologies, 
SFRs, and Schmidt-Kennicutt relations are indistinguishable; see 
A ppendix [Ajl. We make on e small modification to the prescription 
in [Robertson & Kravtsov| ( |2008j l, which is to track the numerical 
pressure floor separately so that it enters into the momentum equa- 
tions, but does not explicitly change the gas temperature (relevant, 
e.g., for determining the cooling function). This is a standard ap- 
proach in high-resolution simulations (see e.g. |Teyssier et al.|20I0| l. 
At the typical resolution we adopt, the pressure provided by equa- 
tion [T] is much less than the turbulent pressure resulting from our 
feedback model (by a factor of ~ 10^ — 10*); only when we turn 
off feedback entirely is the ISM pressure resolution-limited. 

In our simulations, stars are assumed to form from dense gas 
with a constant efficiency e per free-fall time tff = \/37r/32Gp, 
above some minimum threshold po, i.e. 



<SP 3/2 r 

pt — — oc p ' for p > po 
ts 



(2) 



This is numerically implemented by turning gas particles into 
stars stochastically following the calculated SFR (probability p = 
1 — exp{—pt,dt / p), where dt < 100— lOOOyr is the simulation 
timestep and also represents how frequently these values are up- 
dated). Because we wish to resolve the dense regions of star forma- 
tion, we typically set no = lOOcm"^, characteristic of large GMCs. 
The efficiency e is empirically measured in dense star-forming re- 
gions to be « 1 — 2%, roughly constant over a wide range of densi- 
ties ~ 10 — 10 cm^ ( jKrumholz & Tan|2007| ; we adopt a canon- 
ical value of e = 1.5% (see also |Leroy et al.||2008[ l. We discuss 
variations about these fiducial choices in ^4.1 [ below. 



2.1 Cooling and Star Formation 

In order to resolve the formation of very dense clumps, we extend 
the standard atomic-l-metal line cooling curves in GADGET (which 
cut off when the gas becomes neutral at < IO^'K) to allow cooling 



^ Recalibrating our "baseline" A{T) at n = lOOcm"^ gives indistinguish- 
able results. The difference (modulo the standard dependence) is much 
smaller than more dramatic cooling curve variations we consider among 
other numerical tests in § |A] which all produce nearly identical results be- 
cause, in all cases, the cooling time is much less than the dynamical time. 
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2.2 Stellar Feedback 

For the reasons summarized in ^ we model stellar feedback by 
depositing momentum into the gas around young star clusters. This 
in turn drives strong turbulence in the ISM. For standard IMFs (e.g., 
|Kroupa|2002| l, the momentum supplied to the ISM by stellar winds, 
supernovae, and the luminosity of young stars are all comparable 
( [Leitherer et al.|1999[purray et al.|2005l >. If supernovae undergo 
a significant Sedov-Taylor phase, the P-dV work done can increase 
their momentum by a factor of ~ 10 (e.g., Thornt on et al.|1998] l. 
Likewise, if the ISM is optically thick to the infrared radiation pro- 
duced when dust reradiates stellar UV photons, the radiation energy 
density builds up, increasing the radiation pressure force by a fac- 
tor of the infrared optical depth t/r. Modeling these processes in 
detail is a daunting task and one that is beyond the scope of this pa- 
per. Instead, we explore the general properties of models in which 
turbulence driven by momentum-deposition is the dominant stel- 
lar feedback mechanism. This is a plausible approximation even 
in the Milky Way, since the hot ISM contributes only ^ 10% to 
the ISM pressure ( [Boulares & Cox|1990[ l. Moreover, in the well- 
studied star forming region 30 Doradus in the LMC, observations 
directly implicate radiation pressure as the dominant mechanism of 
stellar feedback ( [Lopez et al.|201l"| . We stress, however, that this 
is not intended to represent a complete model of stellar feedback 
and the ISM; in future work, we will study how galaxy properties 
are further modified with the addition of other mechanisms such as 
SNe and stellar wind shock-heating and mass loss, and photoion- 
ization heating. 

In order to make our simplified feedback model as realistic as 
possible, we implement the feedback so that it is explicitly associ- 
ated with young star clusters. We do so by identifying star-forming 
clumps and depositing momentum into the surrounding gas radi- 
ally away from the star clusters. In the following subsections we 
describe the key steps in this method. 



2.2.1 Star- Forming Chimps: Identification 

The first step is to identify star forming regions or nascent star clus- 
ters in the simulation. Starting from each gas particle, we identify 
the nearest dense gas "clump" by iteratively performing a friends- 
of-friends search with an adaptive linking length. Specifically, we 
search over all gas particles within a radius A'smi /Jsmi (with typical 
A'smi ~ 3) of the initial particle to find that particle with the high- 
est local density, and iterate either until a higher-density neighbor 
is not found or until some maximum cutoff is reached. For the lat- 
ter we adopt a maximum of 20 iterations or a distance > 20 times 
the initial particle hsmi (in practice, this limit is rarely reached, but 
is necessary to prevent cases where, e.g., the linking chain might 
traverse a large fraction of the length of a spiral arm). Some care 
is needed in choosing the appropriate value of Mmi, based on the 
physical scales that are or are not resolved in a given simulation ~ 
for our resolutions, A'smi < 1 will simply return the local gas parti- 
cle, and A'smi S> 5 tends to over-link clumps in dense regions such as 
spiral arms and galactic nuclei. Our experiments show that within 
the range A'smi ~ 1.5 — 4 the identification of the peak local density 
is converged for > 90% of all "clumps" (with the remainder mak- 
ing little difference in global quantities, as we show explicitly in the 
Appendix); the density peaks identified in this way also agree well 
with visual identification of overdensities. We thus adopt a canoni- 
cal value of A'smi ~ 3. 

This friends-of-friends search defines the star-forming clump 



of which the initial gas particle is a member.'' The distance between 
the original particle and the clump center (the clump density peak) 
defines the "clump radius" Rdump (if this distance is less than twice 
the initial smoothing length, we set it to this minimum value, since 
the "clump" is effectively unresolved). The enclosed "clump mass" 
in gas (Mciump.gas) or stars (Mdump.*) are defined as the mass of each 
component within a distance of iJciump of the clump center. 



2.2.2 Momentum-Loading 

In our model, stellar feedback is tied to the properties of the stars 
in the stellar cluster associated with a given gas particle. Moreover, 
we only apply the feedback to gas particles that are within 3 /ismi of 
a young star particle (typically < lOpc). This helps ensure that the 
feedback is spatially correlated with young stars. ^ We now motivate 
our implementation in terms of feedback by radiation pressure on 
dust grains. 

At each timestep, we identify the stars (of those formed since 
the beginning of the simulation) within the previously-identified 
clump, and sum their bolometric luminosity, which is a function of 
the star's (known) stellar age 



(3) 



We tabulate L^/M, as a function of age using a Starburst99 iLei- 
[therer et al.|1999^ single stellar population with a [Kroupal p002| > 
IMF at solar metallicity (this time dependence can be important on 
GMC timescales, in contrast to models where all energy is coupled 
instantaneously; see e.g. [Slyz et al.|2005^ . Given the uncertainties 
in the mass-loading factors below, and the fact that our initial con- 
ditions are all relatively evolved systems, it makes little difference 
whether we explicitly allow for a metallicity dependence.* Assum- 
ing that the stellar flux is equally distributed among all of the gas 
within i?ciump, we obtain the luminosity Lj incident on the particle 
in question, which has a mass Mgas,;: 



Lj — Ltot(< ^clump) 



as. y 
(< ^clump) 



(4) 



Because the luminosity incident on a particle in this simple formu- 
lation depends on the light-to-mass ratio of the surrounding stars, 
we find that our results are relatively insensitive to whether we use 
the starlight within i^ciump or some multiple of this radius. 



* We have also experimented with using the center of stellar light or the 
stellar density peak as the location of the clump center (see the Appendix 
for details). In the systems we model here, there is no detectable difference 
between these choices and our fiducial choice of centering the clump on 
the peak gas density. However, the distinction between peak gas and steflar 
quantities could be more important in systems where the main sequence 
fifetime exceeds the dynamical time (e.g.. galactic nuclei) and massive stars 
may wander away from their nataf GMCs. 

^ Because the momentum deposition falls off for gas further from the stars, 
formally extending this to all of the gas makes no difference to our resufts. 

* We neglect for now the fact that at extremefy high resolution, a < 100 Mq 
"star particle" may not completely sample the stellar IMF, and simply take 
the average L« /M» for the particle age. Since we focus on galaxy-average 
quantities, this is probably not a farge uncertainty. But in low-mass clusters 
and GMCs, a more accurate model - for example one which the stellar mass 
range of each particle is sampled stochastically from the IMF, as discussed 
in |Mashchenko et ar|j2008|| - could give interesting differences. 
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Given the incident luminosity, we take the rate of momentum 
deposition in the gas to be 

Pj = (l+»7„riR)^. (5) 

This is the core equation of our feedback model (with typical val- 
ues of riR in Fig.|5j. This force is directed radially away from the 
clump center (i.e. along the vector Rciump). If the particle j itself 
is the clump center, the direction of the force is randomly chosen 
isotropically. 

The first factor of Lj/c in equation [s] represents the momen- 
tum imparted as the optical-UV photons emitted by massive stars 
are absorbed by dust, which re-radiates in the IR. The factor of 
TmLj/c accounts for the momentum imparted by the total number 
of IR photons absorbed/scattered within the gas parcel. Note that 
riR is the optical depth through the clump, not the optical depth of 
the given gas particle. It is the former that sets the total momentum 
supplied to the gas. Finally, equation |5] includes a dimensionless 
factor rjp ~ 1 that accounts for other sources of momentum and un- 
certainties introduced by our simplified treatment. Note, e.g., that 
we do not explicitly include the momentum deposited by stellar 
winds and supemovae separately from that due to the radiation of 
massive stars; rjp > 1 crudely accounts for these additional contri- 
butions.' On the other hand, rjp < I might be appropriate if pho- 
tons efficiently leak out through holes in the gas distribution (see 
Appendix [B|). 

Why do we associate the feedback with the clump and direct 
it from that center of density, as opposed to simply identifying it 
with each star individually? Recall, we are modeling the effects of 
radiation pressure in the limit in which the gas is at least somewhat 
optically thick. If the UV/optical photons could free-stream, then 
the appropriate sources would indeed be each star particle. How- 
ever, if a number of stars are embedded in a gas clump, then in the 
limit of large optical depth, all of the stellar luminosity is trapped 
and re-radiated, so that the momentum flux is everywhere the full 
cItL/c directed radially from the clump center of density and fol- 
lows the scalings we adopt here. This is trivially true in spherical 
or cylindrical (filamentary) geometries, but is a good approxima- 
tion even for complex density distributions if the optical depths are 
sufficiently large. This is an important distinction that makes radia- 
tion pressure a particularly efficient feedback mechanism in dense 
regions (relative to other sources of energy or momentum such as 
SNe or stellar winds). In Appendix[B]we discuss the more compli- 
cated case of an inhomogeneous density distribution. However, to 
the extent that it modifies our conclusions, it is usually equivalent 
to variations in the net efficiency (encapsulated in rjp), rather than 
the spatial distribution or direction of the flux. Of course, if desired, 
the momentum could be isolated to each star by simply taking the 
limit Mmi 0} 

Given that the local density structure of the gas is at least par- 
tially resolved, we use this information to estimate the IR optical 



' We have considered experiments where we include a separate, explicit pj 
term for the direct momentum flux from stellar winds and SNe ejecta, with 
both tabulated in STARBURST99 as a function of stellar age. The absolute 
magnitude of these is, for a constant SFR, ~ L/c. We find this make no 
difference compared to equivalent variation in rjp. 

^ In Appendix [a] we show that directing the momentum from the cloud 
density peak, center of gas mass, or center of stellar mass or luminosity 
makes no difference to our conclusions. Likewise allowing for more com- 
plex geometries by directing momentum along the local density gradients 
gives nearly identical results. 



depth riR = Eeff kir where Eeff ~ Afdump/7i-/;^,„,„p is the gas surface 
density of the clump of interest.' The opacity at IR wavelengths 
is approximately constant for dust temperatures ~ 100 — 1000 K, 
so we adopt kir ~ 5cm"g~' (this is convenient given that we are 
not performing radiative transfer and thus do not have information 
about the true dust temperature). Note that both the weighting of 
Lj and this calculation of tir implicitly scale so that gas near the 
density center where the flux and optical depth are largest will be 
more strongly accelerated than gas in the outskirts of the system. 

We can apply the force associated with pj from equation |5] 
in two ways, either stochastically or as a continuous acceleration 
(the latter is the simplest to implement in grid-based calculations). 
In the stochastic model, we model the momentum deposition by 
randomly "kicking" particles, with an average mass flux set by 

MwVn-^Pj (6) 

where M„ is the mass-loading and i'„ is the initial velocity. What is 
the appropriate "initial" velocity v,v? Models of momentum-driven 
outflows argue that gas should be accelerated to the local escape 
velocity from the star clusters and/or clouds from which they are 
launched ( [Murray et al.|20l01 >. We therefore take 

Vw ~ Vesc « Vv ^dim (Mcliimp , P, ■■■) (7) 

where Vdim is an estimate of the escape velocity as a function of 
the simulation parameters and rjr is a normalization parameter that 
accounts for details such as the exact mass profile shape, micro- 
physical acceleration as a function of position, etc. In practice, we 
have experimented with a variety of choices for the velocity and 
will show that it makes relatively little difference. This is because 
the key parameter that determines the effect of the feedback is the 
total momentum/force (eq.|5]l. 

The escape velocity from the clump as resolved by our simula- 
tions is Vdim ~ x/2GMciump/i?dump. However, some fraction of the 
clump will turn into stars in a dense stellar cluster, the internal dy- 
namics and peak density of which are unresolved. The true relevant 
escape velocity from the location where the outflows are driven is 
probably the escape velocity from that cluster. We therefore take 
the mass in young stars in the clump to be the "star cluster" mass 
and use the empirical size-mass relation of clusters (e.g., [Murrayl 
|2009[ l to determine the cluster escape velocity: 

forM..ci ~ 10^ - 10' Mq. In our models, we take Vdim to be the 
maximum of either the resolved clump escape velocity or the in- 
ferred star cluster escape velocity; the latter is almost always larger. 
In a timestep At, the probability that the particle of mass Mgasj is 
"kicked" is then given by 

= 1 -exp{-(M„, AO/Mg,,.j}. (9) 

The particle then has a momentum Apj = Mgas jVw added to its 
initial momentum, directed radially away from the clump center. 

In addition to the stochastic acceleration of particles described 
above, we can alternatively accelerate the particles continuously 
rather than with individual "kicks"; in this case the particle is sim- 
ply given a Av ; = pj At /Mgasj every timestep. Which prescription 

' For a log-normal density distribution within a given clump, the effective 
optical depth of the inhomogeneous clump is typically within 30% of the 
mean optical depth (Murray et al. 2010). Thus using the latter to determine 
Pj is sufficiently accurate for our purposes. See also §[3.3| 
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is more physically appropriate depends on whether the outflows 
generated by stellar feedback are being accelerated at reasonably 
large radii (e.g. at the outskirts of clouds), or whether they are 
launched in the dense central regions near the star cluster. We show 
below that the two methods yield similar results. 

Many implementations of stellar feedback in the literature turn 
off the hydrodynamics, pressure forces, cooling, and/or star forma- 
tion for some period of time, often chosen such that a wind escapes 
the galaxy entirely (or until the wind reaches some distance from 
its launching point; see'Spring el & Hemquist|2003a[|Oppenheimer| 
[&Dave 2008. Sales et al. 2010^ In our models, by contrast, there 
is no such modification of the underlying equations. We are able 
to directly model the feedback and the resulting dissolution of star 
clusters for three reasons: first, our high resolution allows us to 
resolve a multi-phase ISM structure into which outflows can prop- 
agate; second, we identify massive star-forming regions and drive 
outflows coherently from them, rather than randomly within those 
regions; and third, because the feedback is momentum-driven, it 
drives strong turbulence even in dense, highly radiative environ- 
ments. In situations where a lower resolution is inevitable (e.g. cos- 
mological simulations), it may be necessary for numerical reasons 
to modify the methods proposed here in order to maintain an effi- 
cient source of stellar feedback. This will be studied in future work. 

The feedback model used in this paper is ultimately defined by 
the two key parameters 77,, and 77,. (eqs.|5]c&|7j. We will discuss the 
consequences of different choices for these parameters below; we 
take rjp ~ 1 and 77,, ~ 1 as our physically-motivated default values. 

2.3 Galaxy Models 

Our goal in this paper is to study the effects of stellar feedback on 
the ISM structure and star formation in galaxies. We do so using 
idealized models of disk galaxies with initial conditions motivated 
by galaxies in both the local and high-redshift Universe. We do not 
attempt to model the cosmological evolution of these disks, and so 
do not include extended gaseous halos or cold flows from which 
they would accrete. Rather, our goal is to study how a given feed- 
back mechanism will change behavior given a specific set of (ob- 
served) disk properties. The methodology for building the initial 
galaxies follows that described in detail in a series of papers (see 



e.g. [Springel et al.||2005[ |Di Matteo et a l |[2005[ [Robertson et al.| 
|2006[|Cox et al.|2006| |Younger et al. 2008]i! The disks each include 
an extended dark matter halo with an NFW profile jNavarro et al.| 
|1996| , a stellar bulge (typically with a |Hemquist|1990 profile), and 
exponential stellar and gaseous disks. In all of the models, the ini- 
tial vertical pressure support for the gas disk is provided by thermal 
pressure. As we describe in ^and|4] however, this thermal energy 
is quickly radiated away and the system approaches a new statisti- 
cal equilibrium with star formation and stellar feedback determin- 
ing many of the properties of the gas disk. 

The simulations are carried out at several different resolutions: 
the "standard" resolution has a total of ~ 3 x 10* particles, with 
~ 10* particles in the gas-l-stellar disk (the initial bulges are small 
and so have fewer particles - thus most of the remaining particles 
are in the dark matter halo). Our "high" resolution simulations use 
10 times as many particles, reaching « lO' particles in the disk. We 
also have at least one "ultra-high" resolution simulation per galaxy 
model with > 10^ particles in the disk (to our knowledge, these 
are the highest-resolution galaxy-scale SPH simulations that have 
been performed to date). The models are generally all run for ~ 20 
dynamical times at Rc (~ 3 orbital times), but they typically con- 
verge to quasi-steady state behavior after just ~ 4 — 5 fdyn- After this 



the evolution is essentially just slow, steady-state gradual gas ex- 
haustion; we have confirmed this in at least one run of each galaxy 
model run for 5 times longer than the "standard" runs. As described 
below, the spatial and mass resolutions in each of the simulations 
depend on the galaxy model. 

We consider four galaxy models, motivated by z ~ 2 high star 
formation rate galaxies (non major mergers), local low-luminosity 
LIRGs, Milky Way like spirals, and SMC-like dwarf galaxies. The 
basic properties of these models are summarized in Table [T| 
Sbc: This simulation is designed to model an intermediate- 
mass, relatively gas-rich star-forming disk in the local universe 
(e.g. a low-luminosity LIRG with Lboi ~ lO'"^" Lq and M ~ 
1 — lOMQyr"'). The galaxy has a total baryonic mass 1.05 x 



Hem- 



10 Mq, with a bulge having a mass Mb = 10 Mq and a 
|quist| (1990) scale-length a = 350pc; a stellar disk with a mass 
Md = 4 X 10' Mq and an exponential scale-length of Vd — 1.3 kpc; 
and an extended gaseous disk with Mg = 5.5 x IO'Mq and an ex- 
ponential scale-length of rg — 2.6kpc. The stellar disk has a sech^ 
vertical profile with a scale height of 130pc; it is initialized with 
a radial dispersion profile so that the local Toomre 2 = 1 at all 
positions. The gas disk is similarly initialized in vertical hydro- 
static equilibrium with Q — I. The initial vertical support of the gas 
disk is via thermal pressure. The dark matter halo has a virial mass 
Afhaio = 1 .5 X 1 0" Mq , concentration c = 11, and a spin parameter 
A — 0.033, chosen to match the typical concentrations and spins 
seen in cosmological simulations jBuUock et al.|[200T] |Vitvitska] 
|et al.|2002) ; this gives a total stellar-to-dark matter mass ratio sim- 
ilar to that inferred for systems of this mass (e.g., by |Moster et al.| 
|2010| l. The disk is, however, baryon-dominated within the central 
~ 5 — lOkpc, and as such may develop spiral and bar instabilities. 

For this galaxy model, our standard resolution has SPH 
smoothing lengths of ~ 5 — lOpc in the central few kpc of the disk. 
Our high resolution simulations have ~ 2 — 5 pc smoothing lengths 
and particle masses of ~ 1000 Mq, while the ultra-high resolution 
simulations have particle masses of 100 Mq and Ipc resolution in 
the bulk of the disk.'" 

High-z: This model is designed to approximate a massive, high- 
redshift, and strongly unstable disk forming stars at a very high 
rate ~ 100 — 400 Mq yr~ ' , typical of massive disks observed at z ~ 
2 - 4 ( Erb et al. 2006 ; Genzel et al. 2008 ; T acconi et al.|2010| l. The 
galaxy has a baryonic mass 1.07X 10"Mq, with M,, = 7x 10' Mq 
(a = 1.2kpc), stellar disk M^ = 3 x IO'^Mq (rj = 1.6kpc), gas 
diskMg = 7 X 10'" Mq (rg = 3.2 kpc), initialized with stellar scale- 
height 320 pc and 2 = 1 in gas and stars. This gives a typical gas 
fraction of ~ 0.5 throughout the stellar and star-forming disk (with 
a larger HI gas reservoir at large radii). The halo has Mhaio = 1 .44 x 
10'^ Mq with c = 3.5 and a virial radius appropriate for that mass at 
z = 2. The system is baryon-dominated out to ^ lOkpc. The spatial 
and mass resolution in these simulations are somewhat larger than 
in the Sbc simulation because of the larger total mass and spatial 
size of the disk; however, the Toomre mass and length-scale are 
also much larger, so this model is in a relative sense actually better 
resolved than the Sbc model. 

MW: This system is initialized to represent a local, relatively 
gas-poor. Milky- Way like disk. The galaxy has a baryonic mass 

The particular choice of gravitational softening is chosen as a compro- 
mise between matching the minimum SPH softening lengths, minimizing 
discreteness effects (see e.g. [Power et al.|2()03^ , and giving a similar max- 
imum resolvable density in each simulation (~ lO'cm"^) that is much 
larger than the mean GMC density but still below densities where processes 
of individual star formation and detailed thermal physics become dominant. 
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Table 1. Galaxy Models 



Model eg m, ^halo <^ V'max 

[pc] [Mq] [Mq] [kms-'] [Mq] [Mq] [kpc] [Mq] [kpc] [pc] [Mq] [kpc] [Myr] 

Sbc 2.5 130 1.5ell 11 86 1.05elO le9 0.35 4e9 1.3 320 5.5e9 2.6 0.36 22 

HiZ 3.5 1700 1.4el2 3.5 230 1.07ell 7e9 1.2 3el0 1.6 130 7el0 3.2 0.49 12 

MW 2.5 220 1.6el2 12 190 7.13elO 1.5el0 1.0 4.73el0 3.0 300 0.9el0 6.0 0.09 31 

SMC 0.7 25 2.0el0 15 46 8.9e8 le7 0.25 1.3e8 0.7 140 7.5e8 2.1 0.56 45 



Parameters describing our galaxy models, used as the initial conditions for all of the simulations: 

(1) Model name: shorthand for models of a high-redshift massive starburst (HiZ); local gas-rich galaxy (Sbc); MW-analogue (MW); and isolated 
SMC-mass dwarf (SMC). (2) tg: gravitational force softening in our highest-resolution simulations (ultra-high-res). "High-res" sims use twice this 
value. "Intermediate-res" four times this value. (3) m,: Gas particle mass in our highest-resolution simulations (ultra-high-res). "High-res" sims use 
eight times this value. "Intermediate-res" fifty times this value. New star particles formed have mass = 0.5 m,, disk/bulge particles Ki nij, and dark 
matter halo particles 5m,. (4) M\^.^\„: halo mass. (5) c: halo concentration. Values lie on the halo mass-concentration relation at each redshift (z = 
for SMC, Sbc, and MW; z = 2 for HiZ). (6) Vmax: halo maximum circular velocity. (7) Mbai^y: total baryonic mass. (8) Mb: bulge mass. (9) a: 
|Hemquist|^1990^ profile scale-length for bulge. (10) Mj: stellar disk mass. (11) rj: stellar disk scale length. (12) h: stellar disk scale-height. (13) Mg: 
gas disk mass. (14) r^: gas disk scale length (gas scale-height determined so that 2=1)- (15) /gas: average gas fraction of the disk inside of the stellar 
Re (Mg[< Rf,}/ (Mg[< Rg\+ Mi[< ReX))- The total gas fraction, including the extended disk, is ~ 50% larger. (15) (tjy„: gas disk dynamical time at the 
half-gas mass radius. 



of 7.13 X 10^" Mq, a bulge with Mt = 1.5 x 10^" Mq, a stellar 
disk with = 4.73 x IO'^Mq {vd = 3.0kpc), and a gas disk 
with Mg = 0.9 X 10"'Mq (r^ = 6.0kpc). The disk gas fraction is 
fg = 0.05 — 0.10 throughout the disk out to ~ 8kpc. The disk is 
initialized with a stellar scale-height 300 pc and Q — I. The halo 
has Mhaio = 1.6 X IO'^Mq, concentration c = 12 and i?vir appro- 
priate for z = 0. Observations suggest that the Milky Way hosts a 
pseudo-bulge or a bar instead of a classical bulge, so we initial- 
ize the bulge with a spherical exponential profile (r,/ — l.Okpc), 
rather than a |Hemquist| ( |1990| profile, but since the bulge mass is 
small this makes little difference to our conclusions. At our ultra- 
high (high) resolution, the force and mass resolution in the gas are 

2pc (5 pc) and 200 Mq (2000Mq). 
DwarfySMC: This model is initialized to be similar to the inferred 
properties of the SMC (before entering the MW halo, at least; see 
[Besla et al.|2010[ and references therein), a typical low-mass, gas- 
rich dwarf. The galaxy has a baryonic mass 8.9 x 10** Mq, with 
a bulge having Mi, — lO' Mq (a = 0.25 kpc), a stellar disk with 
Mrf = 1.3 X 10* Mq (rd = 0.7 kpc), and gas disk with Mg = 7.5 x 
10* Mq (rg — 2.1 kpc). The disk is initialized with stellar scale- 
height 140pc and e = 1 . The halo has Mm„ = 2x10"'Mq,c=15 
and Rvii appropriate for z~0. The system is dark-matter dominated 
at all radii outside of the central few hundred pc. For this model, our 
high resolution simulations have a spatial resolution and particle 
mass of < Ipc and ~ 100 Mq, respectively. 



3 GLOBAL GALAXY PROPERTIES 

The key simulations described in this paper are summarized in Ta- 
bles [T] and [2] Table [T| summarizes the properties of each galaxy 
model. Table |2] summarizes the resolution of each simulation, the 
parameters of the star formation model, and the key feedback pa- 
rameters rjp and r/,, (eqs.|5]&|7j. 

Figure [T] shows face-on and edge-on images of the gas sur- 
face density distribution for simulations of each galaxy model with 
our fiducial parameter choices rjp — rj,, = I. Each image is shown 
in the quasi-steady feedback-regulated phase that sets in after a 
few dynamical times. The overall qualitative evolution is similar 
in all of the simulations with feedback. The gas cools efficiently 
to low temperatures and collapses by gravitational instability at the 



Jeans/Toomre scale. This leads to the formation of dense clumps 
that are the sites of star formation and, in our model, feedback. 
The resolved density contrast between the centers of star-forming 
clumps and the interclump medium is typically ~ 1000 but can 
be as high as ~ lO''. The ISM sustains this clumpy structure as 
long as we evolve our simulations, as gas is blown out of individ- 
ual clumps (by feedback) into the more diffuse ISM before being 
incorporated into new dense clouds. We defer a rigorous analy- 
sis of the lifetimes and evolution of individual clumps for future 
work (in preparation) analyzing the properties of GMCs, where we 
can make rigorous comparisons with observations. But typically, 
we find average lifetimes of individual clouds < 10 Myr or a few 
free-fall times (weakly increasing with mass scale from the SMC 
through HiZ models), giving an integrated fraction ~ 1 — 5% of 
clump mass turned into stars. 

This fragmentation is the natural extension of Jeans-mass 
GMCs in the MW and other nearby galaxies. Indeed, if we wish 
to explicitly resolve these scales, most of the gas mass should be 
in dense sub-clumps at something like the Jeans scale. The primary 
role of feedback is to regulate against runaway collapse and star 
formation in those clouds. Figure[2]illustrates how these morpholo- 
gies depend on the strength of feedback. We consider the HiZ case, 
which is most strongly unstable, at two different extremes (holding 
all details of the model fixed, except feedback strength). First, with 
feedback much stronger than is realistic, rjp = 100. In this case, 
essentially all sub-structure in the galaxy is "wiped out," and the 
star formation is smoothly distributed over a ~ 10 kpc disk (despite 
< lOpc resolution). This would be analogous to a MW model with 
no GMCs, where all star formation occurred in regions with local 
densities ~ 1 cm^'. Second, we consider a case with no feedback. 
In this extreme, the opposite occurs: the GMC complexes seen in 
Figure |2] dissipate their internal velocity dispersions and experi- 
ence runaway collapse and star formation. This collapse proceeds 
until the GMCs reach the simulation resolution limit and leads to 
all of the gas being at extremely high densities, n ~ 10* cm""' (as 
we show explicitly below); a corollary is that the gas is converted 
into stars on essentially one (large-scale) dynamical time. 

Figure [3] shows the phase diagram for the gas in each of our 
fiducial simulations: we plot both the thermal sound speeds and 
turbulent velocities as a function of gas density (averaged over the 
smallest available scale, the SPH smoothing length). For low den- 
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Figure 1. Images of the gas distribution for our fiducial simulations {rjp = r)v = 1) in the feedback-regulated quasi steady-state. Brightness shows the gas 
surface density while color shows the specific SFR (increasing from blue to red); both are on a logarithmic scale spanning a dynamic range of ~ 10*. 
Top: Large scales (wide-field image) out to twice the half-gas mass radius. Middle: Intermediate scales (zoom-in of the image at top) out to the half-SFR 
radius. Bottom: Edge-on; scale is same as the middle image. One example is shown for each of the initial conditions we model (HiZ_10_4_hr, Sbc_10_4_hr, 
MW_10_4_hr, and SMC_10_4_hr in Table|2}. The simulations develop complex substructure and exhibit a diverse range of gas morphologies. Most stars are 
formed in dense but resolved giant 'molecular' cloud complexes, which are the sites of the feedback modeled here. 



HiZ: Extreme Feedback ni^.: no reeocacK 

^ ■ ■ ■ 

1 kpc ^^^^ 1 kpc 



Hi?: No Feedback 



idback 



Figure 2. As Figure^(m(drf/e left), but for an otherwise identical HiZ sim- 
ulation with extremely strong feedback {left) with rjp = 100 (this is not a 
realistic choice but purely shown for illustrative purposes), and one with no 
feedback (right). With arbitrarily strong feedback, all collapse of gas into 
GMC complexes is suppressed. With no feedback, the cloud complexes in 
Figure^undergo runaway collapse to the resolution limit (the single white 
pixels at right); the mass piles up at densities > lO* times larger than in our 
"standard" models. 



sity gas with n <C lcm^% the sound speed and turbulent velocity 
are often comparable, but for denser gas the turbulent velocity is 
always much larger than the thermal sound speed. 

The characteristic densities of clumps/GMCs are evident in 
the peak of the gas distributions near n ~ lOOcm"^ in Figure [s] 
the typical turbulent mach numbers for this gas are ~ 30 — 100. 
Because of the high Mach numbers, turbulent motions rather than 
thermal motions are the dominant impediment to gravitational col- 



lapse. Specifically, the characteristic mass of large GMCs is set by 
the turbulent Jeans mass for the bulk of the matter, and corresponds 
to: ~ 10^ M0 in the SMC case, ~ i(f Mq in the MW and Sbc 
cases, and ~ 10* Mq in the HiZ case. These estimates agree rea- 
sonably well with the observed properties of massive cloud com- 
plexes in the respective systems. By contrast, if the gas were ther- 
mally supported, the characteristic mass of collapsed gas would be 
much smaller. For the dense gas, however, thermal support is only 
important on scales below the sonic length (< 0. 1 pc), which is well 
below our resolution limit. 

The minimum pressure to prevent unresolved collapse below 
the resolution limit (eq.[TJ is well below the resolved turbulent pres- 
sure for the median densities in Figure [3] This effective pressure 
does, however, produce the small "upturn" in the turbulent 5v at 
the very highest densities n ^ 10"* cm^'. For our purposes, the key 
point is that we resolve the median GMC length, density, and mass 
scales well, even in our lowest-resolution models. 



3.1 Morphologies 

There are a variety of morphologies present in the simulated galax- 
ies depending on how self-gravitating the disk is. The high-redshift 
disk analogues (HiZ) are the most strongly self-gravitating and 
so fragment into very massive clumps (Mxoomrc ~ 10* — lO' Mq), 
which dominate the star formation. This morphology resembles the 
clumpy systems observed at z ~ 2 — 3 jGenzel et al.|2008{|Taccoml 
|eraLl[2006| [Law et a l. '2009'). This is even more clear when we 
focus on the region which contains half the star formation {mid- 
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Table 2. Simulations Plotted in This Paper 



Simulation 




SF Law 






HiZ_8_0_nofb 


2e6 








HiZ_8_2_nofb 


2e7 




- 


— 


HiZ_10_4 


2e6 


_ 


1.0 


1.0 


HiZ_8_10 


2e6 




2.0 


2.0 


HiZ_8_ll 


2e6 




4.0 


2.0 


HiZ_9_l 


2e6 




10.0 


I.O 


HiZ_6_0_hr 


le7 




1.0 


1.0 


HiZ_7_0_hr 


6e7 


_ 


1.0 


1.0 


HiZ_6_3_hr 


le7 


p oc p'" 


1.0 


1.0 


HiZ_6_4_hr 


le7 




1.0 


1.0 


HiZ_7_l_hr 


6e7 


e = 0.35% 


1.0 


1.0 


HiZ_7_2_hr 


6e7 


e = 6.0% 


1.0 


1.0 


HiZ_7_3_hr 


6e7 


«c = 2500 


1.0 


1.0 


HiZ_7_4_hr 


6e7 


fJc = 25 


1.0 


1.0 


HiZ_10_4_hr 


2e7 




1.0 


1.0 


HiZ_10_5_hr 


2e7 




1.0 


2.0 


HiZ_10_6_hr 


2e7 


- 


1.0 


4.0 


HiZ_10_7_hr 


2e7 


- 


2.0 


1.0 


HiZ_10_8_hr 


2e7 


- 


4.0 


1.0 


HiZ_10_9_hr 


2e7 


- 


10.0 


1.0 


HiZ_10_14_hr 


2e7 


- 


0.33 


1.0 


HiZ_10_ll_hr 


2e7 


- 


1.0 


a 


HiZ_8_14_hr 


2e7 


- 


4.0 


2.0 


HiZ_8_17_hr 


2e7 


- 


5.0 


4.0 


HiZ_10_4_uhr 


2e8 


- 


1.0 


1.0 


MW_8_3_nofb 


3e6 








MW_9_1 


2e6 




1.0 


1.0 


MW_10_7_hr 


le7 




1.0 


1.0 


MW_10_8_hr 


le7 


e = 0.35% 


1.0 


1.0 


MW_10_9_hr 


le7 


e = 6.0% 


1.0 


1.0 


MW_10_10_hr 


le7 


poc p" 


1.0 


1.0 


MW_10_ll_hr 


le7 




p oc p- 


1.0 


1.0 


MW_10_12_hr 


le7 


«c = 10 


1.0 


1.0 


MW_10_13_hr 


le7 


= 1000 


1.0 


1.0 


MW_9_l_hr 


2e7 




1.0 


1.0 


MW_9_2_hr 


2e7 




1.0 


a 


MW_9_3_hr 


2e7 


- 


1.0 


2.0 


MW_9_4_hr 


2e7 


- 


0.33 


1.0 


MW_9_5_hr 


2e7 


- 


10.0 


1.0 


MW_8_4_hr 


3e7 




10.0 


4.0 


MW_8_5_hr 


3e7 


- 


10.0 


1.0 


MW_10_2_hr 


3e7 


- 


10.0 


2.0 


MW_10_4_hr 


3e7 




1.0 


1.0 


MW_9_l_uhr 


2e8 




1.0 


1.0 


MW_8_uhr 


3e8 




10.0 


2.0 


SMC_10_3_nofb 


2e7 


- 


- 


- 


SMC_10_l_hr 


2e7 


- 


4.0 


2.0 


SMC_10_2_hr 


2e7 




10.0 


2.0 


SMC_10_4_hr 


2e7 




1.0 


1.0 


SMC_10_uhr 


le9 




10.0 


2.0 


Sbc_10_3_nofb 


2e7 








Sbc_10_l_hr 


2e7 




4.0 


2.0 


Sbc_10_2_hr 


2e7 




10.0 


2.0 


Sbc_10_4_hr 


2e7 




1.0 


1.0 


Sbc_10_uhr 


2e8 




10.0 


2.0 



Parameters of our key simulations (only simulations appearing 
in Figures are listed; others are noted in the text): 

(1) Name/ID. First characters correspond to the class of galaxy 
model ("SMC," "MW," "Sbc," or "HiZ," as in Table[T). 

(2) Total particle number 

(3) Star formation law. "-" corresponds to the default law: 

p* = ep/%(p) for p > po, with e = 1.5% and«o = 100 cm~^; 
varied quantities are noted (see Fig.[6|. 

(4) Momentum-loading normalization rjp (see eq.|5j. 

(5) Initial velocity normalization r;,. (see eq |7|- 

© cfcdffKte?*9NKASn(}0B}<fHK)Ci0ter than discrete "kicks." 



100 



10 



HiZ 




MW 



Turbulent o(hj.^|) 
Thermal 





SMC : 



2 4 
Log(n) [cm"] 



6 -2 2 4 6 
Log(n) [cm "] 



Figure 3. Phase diagram for the gas in the fiducial simulations in Figure^ 
at times in the feedback regulated quasi steady-state. Contours are iso- 
density at ~ I0~'' -2,-1.5.-1, -0.5 maximum (progressively darker 
dotted, short-dash, dot-dash, dot-dot-dot-dash, long-dash, solid contours, 
respectively). Blue contours show the thermal sound speeds c, oc r'/-; 
red contours the local turbulent velocity dispersion a (averaged within one 
gas smoothing length /Zji^i around each particle). Lines of constant Jeans 
mass oc (5v^^h~'/^ (black dotted) are shown for comparison. The median 
clump/cloud gas density is evident in the peak near ~ lOOcm"^. For all 
the dense gas, the thermal pressure is negligible compared to the turbu- 
lent pressure/velocities. As a result, the turbulent Jeans mass governs large- 
scale collapse and corresponds to the mass of massive clumps/GMCs (from 
~ 10^ Mq in the SMC model through ~ 10** Mq in the HiZ model); these 
are very well-resolved. The thermal Jeans masses are much smaller, but are 
only relevant for the dynamics on scales below the sonic length (< 0. 1 pc) 
where individual groups of stars form; this is unresolved, hence the neces- 
sity of an "effective" small-scale star formation law. The slight "upturn" in 
o"('isml) at n 2> 10*cm~^ (most evident in the HiZ model) comes from the 
minimum pressure corresponding to the |Truelove et al.|^1997} Jeans condi- 
tion (eq.^. This indicates where resolution limits prevent us from following 
further collapse to higher densities. 



die panel) ~ this is dominated by a few giant complexes. Viewed 
edge-on, the HiZ model appears qualitatively similar to the "clump 
chain/cluster" systems observed at high redshift. 

The Sbc model fragments in a manner similar to that of the 
HiZ model. However, the disk is thinner, and the Jeans mass and 
length-scales are significantly smaller, so star formation is more 
distributed in many clouds (the number of massive clouds predicted 
at 2 ~ 1 is ~ {R/hY and is thus larger for thinner disks). At slightly 
later times than that shown in Figure[T| the system develops a strong 
stellar bar, and the gas - while still very clumpy - flows into the 
center along the m — 2 mode. Flocculent spiral structure also de- 
velops in some of the Sbc runs at large radii. 

The Milky Way model in Figure [T| shows clear grand design 
spiral structure. It is also weakly barred in the center, though the bar 
feature is much more prominent in the stars. Note that this model 
has a lower gas fraction and is significantly more stable than the 
Sbc and HiZ models - the latter because it is dark matter or bulge 
dominated at all radii. As a result, the characteristic clump mass is 
smaller (Mioomrc ~ 10^ Mq) and the Jeans length is much smaller 
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Figure 4. Total star formation rate for each of our galaxy models in Table[T]as a function of time, both with feedback (tj,, = '^r = 1 ) and without. The timescales 
are different in each model and correspond to the characteristic dynamical timescales in each system (longer in the more stable, dark-matter dominated systems; 
see TableQ. Absent feedback (red dot-dashed line) the gas collapses on a dynamical time, leading to a SFR well in excess of that observed in similar systems; 
the SFR then declines as the gas is exhausted. With stellar feedback, the SFR reaches an approximate equilibrium in which feedback maintains marginal 
stability to gravitational collapse (Q ~ 1). 



relative to the effective radius (for a 2 ^ 1 disk, Ajeans ~ fgnsRe, 
so ~ 100 times smaller here). As a result, the individual "clumps" 
are much less prominent in the image, despite the fact that most 
of the mass in the star-forming disk does lie in thousands of re- 
solved "clouds" with masses ~ 10"* - 10* Mq. The small gas frac- 
tion also causes the disk to be significantly thinner in the edge-on 
image: Q ~ 1 implies h ~ fgR for weakly self-gravitating disks 
(e.g., [Thompson et al.|2005] l. In future work (in preparation), we 
will investigate the detailed structural properties of the ISM and 
simulated GMC-analogues to compare them to observations of the 
MW and local group galaxies. 

The SMC-like model behaves quite differently from the Milky 
Way model, although both are dark-matter dominated. The SMC 
model is completely stable to global instabilities and thus forms 
stars in a more uniformly distributed fashion. The ISM on these 
scales is turbulent and patchy, with an irregular or (on large scales) 
featureless structure, typical of observed dwarf galaxies. Despite 
the low SFR of ~ O.lMQyr"', the turbulent velocities generated 
by stellar feedback are sufficient to make the system quite "puffy" 
and thick (given the weaker potential depth). Figure [T] shows that 
individual star-forming regions are resolved with size scales of < 
lOpc. 

Note that because the gas in this model is quite low-density, 
the cooling times are long and energy input via supemovae and 
stellar winds will have a significant effect on the gas morphology. 
There are plain indications here that the present model, including 
momentum from radiation pressure alone, is not a complete de- 
scription of the ISM. For example, the temperature of the "diffuse" 
ISM in all the galaxy models tends to be much too low. We show 
this explicitly in Figure [3] where we plot the phase distribution of 
the gas. The volume-filling gas distributed between dense clouds 
is almost entirely "warm" (lO'* < lO' K), with negligible mass 



in the characteristic "hot phase" of the ISM atT> 10* K (there is 
some, generated by shocks, in the stronger HiZ and Sbc cases, but 
even here it is less than a percent of the total gas mass). Some ad- 
ditional heating mechanisms, such as SNe and "fast" stellar winds, 
are probably critical to explain the full temperature structure of the 
ISM. In future work we will investigate this in detail, with explicit 
models for various heating terms; for now, we simply note that the 
small mass fraction in the "hot" phase, while potentially important 
for phenomena such as galactic winds, is unlikely to change the 
structure of cold regions as it contains little mass and, even in MW- 
like galaxies, contributes only ~ 10% to the typical ISM pressure 
( [Boulares & Cox|1990| l. We see in Figure[3]that the turbulent veloc- 
ities are much larger in all dense gas than the thermal sound speeds 
(and tend to be near-virial), making the detailed thermal structure 
sub-dominant on these scales. 

3.2 Star Formation Histories 

Figure |4] shows the star formation history (galaxy-integrated star 
formation rate [SFR] as a function of time) of each of our galaxy 
models for the same feedback parameters used in Figure[T] we also 
compare to simulations of the same galaxy models that include 
cooling and star formation, but not stellar feedback. 

In the models without feedback, the SFR increases to a peak 
value on a single global dynamical time; the SFR remains at this 
value until the gas in the disk is exhausted. The peak SFRs in the 
simulations without feedback are a factor of > 10 larger than those 
observed in the systems that motivate these galaxy models - the 
observed values are ~ (50 - 300, 3 - 20, 2 - 4, 0. 1 - 0.5) yr" ' 
for high-z non-merging SMGs ([Forster Schreiber et al.|2009^, low- z 
non-merging LIRGs (Sanders &' Mirabel|19 96"Evans et al.|2009a| l, 
the MW and similar-mass spirals at z = 0, and isolated SMC-mass 
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systems at z = jNoeske et al.|2007[[Salim et al.|2007| >. In ^we 
explicitly show that these models also lie well off of the observed 
Kennicutt-Schmidt relation between SFR and gas surface density. 
Physically, this is because in all of our simulations the gas can cool 
to an arbitrarily low temperature on a timescale short compared to 
the local dynamical time. In the absence of feedback, the gas is 
then violently unstable to runaway clumping and rapid star forma- 
tion. The net result is that M, ^ Mgas/?dyii, i.e., most of the gas is 
converted into stars on a single dynamical time, the timescale for 
the initially thermally supported gas disk to collapse. This behavior 
is physically correct in the absence of stellar feedback, and should 
be recovered in any simulation that does not include such feedback. 

A number of authors have suggested that instabilities due to 
self-gravity alone might generate the turbulence needed to slow 
down star formation in g alaxies (e.g. [B allesteros-Paredes et al. 
[2007{ [Tasker & Tanl[2009l |Krumholz & B urkert 2010 1. Figure |4 
is not consistent with this hypothesis. Absent stellar feedback, the 
majority of the gas accumulates into dense clumps in which star 
formation proceeds unimpeded. Independent simulations at simi- 
lar resolution but with different physics included have reached the 
same conclusion (e.g. |Boumaud et al.|2010| ). We thus find that stel- 
lar feedback is critical to regulating star formation in galaxies. A 
more subtle question is, when strong feedback is present, does it 
"drive" the turbulence, or is it still primarily driven by gravity? We 
will investigate this more quantitatively in future work. It gener- 
ally appears, however, that the role of feedback is to offset the dis- 
sipation of turbulence and relative motions (particularly in dense 
regions), so in this sense it "provides" momentum; but the level 
it must provide, and the turbulent cascade and regulation of those 
motions, is primarily dictated by gravity. 

In contrast to the models without feedback, our simulations 
with stellar feedback rapidly reach a maximum SFR and then 
remain at this quasi-steady state for many dynamical times. In 
some of our simulations, there is a slow increase in the SFR on 
a timescale longer than the disk's dynamical time; this is a conse- 
quence of both slowly-growing secular instabilities (e.g. halo bars) 
and spatial redistribution of gas by the stellar feedback (e.g. gas 
being driven in fountains from small to large radii). Eventually, 
however, since these are isolated systems without continuous gas 
accretion, the SFR must decline by gas exhaustion, but this decline 
is much more gradual than in the absence of feedback. Test runs, 
run for ~ 5 times longer, confirm that there is no new behavior after 
the first few dynamical times; the SFRs gradually decline as gas is 
exhausted. In ^we show that our simulations with stellar feedback 
are reasonably consistent with the the observed Kennicutt-Schmidt 
relation. This is true for a range of feedback parameters. Below we 
discuss the physical origin of this low star formation efficiency. 

3.3 Structural Properties 

Figure [5] shows a number of the properties of the ISM in our HiZ 
model as a function of time and radius: the vertical velocity dis- 
persion {a-: thermal and turbulent), the Toomre Q parameter" of 
the gas, the mass weighted density distribution function, the total 
momentum supplied to the gas, the "kick" velocity particles ini- 
tially receive, and the (momentum-weighted) optical depth of gas 
clumps where the kicks are applied (i.e., the optical depth for the 

" We define Q here as crK/TrGSgas, where a is the full gas velocity dis- 
persion, K is measured from the azimuthally-averaged mass profile, and 
Sgas is the gas surface density. 



regions where most of the feedback occurs). These results provide a 
more quantitative view of the quasi-steady feedback-regulated state 
reached in our simulations. We show the results for the HiZ model 
because the strong gravitational instability, high gas fraction, and 
very high star formation rate make it the model most sensitive to 
variations in our feedback prescription and give the largest differ- 
ences between models with and without feedback. However, we 
find identical qualitative conclusions (discussed below), modulo 
the absolute value of the various quantities, for each of our other 
galaxy models. We focus on three different simulations in Figure|5] 
The first is one of our ultra-high resolution runs (HiZ_10_4_uhr) 
with r),, = 77,. = I and 2 x 10* particles, in which a typical Jeans- 
mass clump in the disk is resolved with as many as ~ 10^ particles. 
We compare this to a lower resolution simulation with the same 
feedback parameters (HiZ_10_4) and to a lower resolution simu- 
lation which has rjp = 10 (HiZ_9_l) to compensate for the poorer 
resolution of the densest star-forming regions. 

Perhaps the most important result in Figure |5] is that the ISM 
properties do not depend sensitively on either resolution or the mo- 
mentum feedback parameter rjp (the star formation history does de- 
pend mildly on rjp as we show in ^4.2[ >. The key reason for this is 
that the disk always self-regulates to maintain 



where 5v is the turbulent velocity dispersion induced by the stel- 
lar feedback. Figure |5] shows explicitly that all of the simulations 
maintain Q ~ 1 in the feedback-regulated phase (top middle panel). 
The differences between models are small and all within the range 
of random variations and noise. Figure |5] also shows the (mass- 
weighted) vertical velocity dispersion a- = ^ c? + 5\'l of the gas 
as a function of time {top left panel). Initially decreases rapidly 
as the thermal support is radiated away. As star formation com- 
mences, however, stellar feedback quickly drives the turbulent ve- 
locity to Sv, ~ 30 — 50kms^'. Given this turbulent velocity, the 
vertical scale-height of the disk is a few hundred pc, with only a 
modest dependence on radius; at all radii this thickness is much 
larger than the resolution limit. 

The early-time and late-time values of a, in Figure|5]are com- 
parable because in both limits 2 ~ 1. The models are initialized 
with thermal support and 2=1 but this is quickly replaced by 
turbulent support that self-consistently maintains 2 --^ 1 at later 
times. The velocity dispersions in Figure [5] are also in reasonable 
agreement with the observed values in high-redshift disks ( ,Forster| 
[Schreiber et al. [2006) 1. The other galaxy models also self-regulate at 
2 ~ 1. However, given their lower masses, gas fractions, and star 
formation rates, this translates to lower absolute velocity disper- 
sions: Sv ~ lOkms^' in the MW and Sbc models, and ~ 6kms~' 
in the SMC model (modulo rescaling by this absolute value, how- 
ever, the dependence of ct, on time, resolution, and rjp is nearly 
identical to that shown in Figure[5](. 

The top right panel in Figure |5] shows the mass-weighted gas 
density distribution averaged over the entire galaxy once the star 
formation reaches an approximate steady-state (since most of the 
gas mass is near ~ 3 kpc, the density distribution in an annulus at 
this radius is quite similar); the distribution is close to lognormal in 
all of the simulations with a median density of ~ lOOcm^^ and a 
broad dispersion of ~ 1.5 dex. The highest resolved densities reach 
> 10*cm~' in the ultra-high resolution simulation, but it is im- 
portant to note that gas does not simply "pile up" gas at these high 
densities, which it does if we do not include feedback. We show this 
explicitly in Figure|5]by including the density PDF for a simulation 
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Figure 5. Properties of the ISM and feedback in several of our HiZ simulations: intermediate (HiZ_10_4) and ultra-high resolution (HiZ_10_4_uhi') simulations 
with Tjp = \ and rjy = 1, and an intermediate-resolution (HiZ_9_l) simulation with ij,, = 10 (see Table|2j. Top Left: Vertical gas velocity dispersion, a- = 
y'c^ + Sv^ (averaged over the entire disk, weighted by gas mass). The initial disk is thermally supported, but this thermally energy is rapidly radiated away; at 
later times a comparable is produced by feedback-driven turbulence. Top Center: Gas Toomre Q parameter in narrow radial annuli as a function of radius 
(averaged over times > 60Myr, when the system is quasi-steady state). Top Right: Gas density distribution (gas mass per interval in log n) is lognormal with 
~ 1 — 1.5 dex dispersion; low and high-res simulations converge to the same median density, but at low resolution the full width is not resolved. With no 
feedback (dotted), the gas piles up at the highest resolvable densities. Bottom Left: Sum of all momentum (| Ap|) injected via feedback (solid; eq.jsj compared 
with input optical-UV stellar photon momentum = J L» c~ ' (dotted). Note that the momentum injected is nearly the same for all three simulations, including 
rip = \ and rjp = 10. The dot-dashed line shows that the input momentum is well-reproduced using the optical depths from the bottom right panel and only 
the very young stars (< 10^ yr old). This demonstrates that star-forming clusters disrupt rapidly. Bottom Center: Mean "kick" velocity given to gas particles at 
their launching from young stellar clusters (and l-cr dispersion); values approach ~ 150— 200kms^' , as expected given the massive 10** Mq clumps forming 
in these simulations. The kick velocity is much larger than the actual dispersion in the disk because the particles shock and share their momentum immediately. 
Bottom Right: Resolved IR optical depths of gas clumps used in the feedback model (eq.|5j. In the simulations with rjp = 1, r ~ 30 — 50, corresponding to 
S ~ 10 g cm^^, comparable to the observed surface densities of star clusters on ~pc scales. The siinulation with rjp = 10 has the same total input momentum 
(bottom left panel) but as a result the gas clumps only collapse to r ~ 10. A comparison of our MW-like models gives identical qualitative conclusions, but 
with systematically shifted absolute values: ^ lOkms"', Q 1, («) ~ 1 cm' -3, "kick" i'~ 30-50kms-', and (r) ~ 10- 30 at r?p = 1. 



with identical initial conditions, but no feedback (see also the den- 
sity distributions in the simulations without momentum-feedback 
in |Teyssier et al.|[2010[ l - in this case almost all the gas ends up 
at the maximum density allowed by our resolution (~ lO'cm^'), 
with a small tail at low densities. With feedback included, how- 
ever, most of the mass is in GMC-like structures, but within those 
structures feedback ensures that most of the mass is in a more dif- 
fuse phase, rather than in the densest star-forming cores. The same 
conclusions pertain to our other galaxy models, but with lower me- 
dian densities as expected; the volume-averaged (n) ~ 1 cm^' in 
the MW and Sbc models, with much of the mass in the star-forming 
disk in GMCs with a mean (n) ~ 10 — 30cm~^ (and a resolved tail 
up to ^ 10*cm~^). We caution that the distribution of low-density 



gas (n <^ 1 cm~') can be strongly altered by other sources of en- 
ergetic feedback, such as SNe, stellar winds, and photo-ionization; 
the most dense gas, however, is where radiation pressure is likely 
to be most important. 

The bottom panels in Figure [5] quantify the magnitude of the 
stellar feedback: we show the integrated momentum supplied to 
the gas as a function of time, the typical initial velocity of the kicks 
as they are imparted to particles at the star-cluster scale, and the 
momentum-weighted optical depth of the gas clumps where the 
feedback is applied. 

The values of the "initial" kick velocities given to the particles 
are large for the HiZ model, but not surprising given the very large 
star cluster masses (~ 10^ Mq) associated with the giant clumps 
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Figure 6. Star formation rate as a function of time for the HiZ model for variations in the small-scale (high-density) star formation law; all runs use our fiducial 
feedback parameters {rip = rjy = 1). These results demonstrate that the global star formation rate depends only weakly on the small-scale star formation law. 
In each panel, the black solid line shows the standard star formation model: p, = ep/tg above a threshold density iiq = 100 cm~^, with e = 0.015 and 
tff = ^^3n/32G p Qc p^"'. Left: Variations in the star formation efficiency t. Middle: Variations in the density-power law of the star formation model: 
p* oc p" with 11 = 1 , 1.5,2, normalized so that p* is the same as the default model at ng. Right: Variations in the threshold density for star formation riQ. 




Figure 7. As Figure |6] but for the MW model. Again the global SFR is 
independent of the local, high-density SF law. 



in gas-rich high redshift systems (see [Murray et al.|2010| [Genzel| 
|et al.|20lT| l; by contrast, the initial kicks in the MW-Uke system are 
much lower, tens of km s ~ ' . Note also that the initial kick velocities 
are much larger than the velocity dispersion in the galaxy (in both 
cases): this occurs because the particles immediately interact with 
the surrounding ISM and share their momentum. In j4.2| we show 
that for this reason, the choice of the initial kick velocity - or even 
whether to continuously accelerate rather than "kick" particles - 
is largely irrelevant.'" Instead, the important parameter is the total 
momentum supplied to the gas. 

With the large SFR and reasonably large kick velocities of the 
HiZ model, we might expect a sizeable super-wind to be generated 
by this feedback mechanism alone. However, in fact, the amount of 
mass in a proper super-wind (defined as e.g. the mass flux at > Vc 
escaping to at least ~ 20kpc) is relatively small relative to the SFR, 
about ~ 10%. This is because, as described above, the momentum 
is coupled in extremely dense regions and so rapidly shared among 
the gas particles. This can maintain a large velocity dispersion in 
the disk, but will not efficiently launch gas well out of the disk. 
Occasionally some material has an un-obstructed sightline out of 
the disk and escapes, but even then, the launch velocities are typi- 
cally below the circular velocity, so the material is lofted up above 

For the same reason, the total fraction of gas particles initially "kicked," 
which in these models is about ~ 10%, is unimportant. 



the disk and then returns rapidly. It is likely that rather than winds 
being launched directly out of the galaxy from individual star clus- 
ters, some continuous acceleration mechanism is needed to act on 
gas once this local mechanism pushes it above the disk, in order to 
accelerate it out of the galaxy halo. This could be either pressure 
acceleration from hot, SNe-heated gas, or continuous radiation ac- 
celeration from the light which escapes the dense, optically thick 
regions we model here. In future work, we will investigate the prop- 
erties of the galactic super-winds in more detail, and examine how 
these mechanisms interact with the feedback mechanism described 
here. 

Figure |5] shows that the total momentum supplied to the gas 
is significantly larger than J (L/c) At = Enijc because of the non- 
zero optical depths (where E^a is the total radiated energy). In fact, 
the numerical results are consistent with the total momentum sup- 
plied being given by ~ rj,, {TiR)£',ad.young/c, where Erad. young is the 
integrated luminosity from young stars with ages < 10* yr. This is 
because the feedback begins to disperse the densest regions on a 
~ lO' yr timescale. Figure [s] also shows the median and disper- 
sion in the clump optical depths for the regions where the feedback 
is applied: tir ~ 50 in the highest resolution simulation. This cor- 
responds to gas surface densities ~ 10 g cm~^, similar to the ob- 
served surface densities of massive star clusters. The average r over 
the entire disk is, of course, significantly smaller, tir ~ 0.1 — 1. For 
this reason, for the MW, Sbc, and SMC models, although the disk- 
averaged TiR is significantly smaller than the HiZ model, their "ef- 
fective" riR is not too much smaller. Despite the global gas mass be- 
ing lower, the mass that actually forms stars and star clusters tends 
to be compact cores at high three-dimensional (n > lO^cm"^) and 
surface densities (E » IOOOMqpc"^). 

Figure|5]shows that the total momentum input does not depend 
that strongly on resolution or on r\p at a given resolution. In fact, the 
optical depths decrease with increasing r\p, maintaining approxi- 
mately the same total momentum input. Physically, this is because 
for more/less efficient feedback the gas collapses to lower/higher 
densities (on average). The fact that the total momentum supplied 
by feedback depends only weakly on resolution and r\p is a con- 
sequence of the disk self-regulating to maintain Q~ \. This con- 
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straint picks out a 5v as a function of Eg and f2 (eq.|10^ - the mo- 
mentum input then adjusts to produce the required 5v. Again, the 
same is true in all galaxy models. 

The column density distribution within individual star form- 
ing clumps can strongly influence the efficacy of radiation pres- 
sure feedback. For example, if there is a very broad distribution 
with a large number of optically thin sightlines, then even though 
the average column density of a clump in the IR may be large, a 
sizeable fraction of photons would leak out of optically thin sight- 
lines and the effective optical depth for the purposes of feedback 
would be reduced {r]p < 1, in our parameterization). To quantify 
this, we considered a number of massive clumps in our HiZ sim- 
ulation and determined the optical depth along ~ 1000 sightlines 
evenly spaced in solid angle from the clump center outwards (fol- 
lowing the methodology in |Hopkins et al.|2005^ . The characteristic 
dispersion in optical depth /or a given clump is small, ~ 0.2 dex, 
similar to what has been found in smaller-scale simulations of indi- 
yidual clouds ([Ostriker et aL|2001) and m easured observationally 
I Wong et al.|2008[ [Goodman et al.|2009| ." Note that this disper- 
sion for an individual star-forming clump is much smaller than for 
the galaxy as a whole. A key consequence of the relatively narrow 
column density distribution within clumps is that only a negligible 
fraction of the sightlines are optically thin enough for the photons to 
rapidly leak out; it thus appears reasonable to use the mean clump 
column density when quantifying the feedback produced by the IR 
radiation (as we assume in our fiducial models; see §2.2.2[ l. 

It is also straightforward to calculate the total momentum de- 
position or infer it from the momentum coupled and typical ve- 
locities, E„ « (1/2) pwV„. For the values in Fig. [s] this is about 
0.2 — 0.4% of the stellar luminosity. Interestingly, recent observa- 
tions of massive stars ( [Freyer et al.|2006) and star forming regions 
( [Lopez et al.|20lf| > have suggested similar ~ 1% efficiencies for 
transfer of luminous energy into bulk motions and (post-shock) 
thermal energy. 

In future work (in preparation), we will compare the structural 
properties of the ISM and dense gas in these simulations and local 
observed galaxies in detail. However, in low-density gas, typical of 
much of the mass in the SMC model and the intermediate/diffuse 
phases in the Sbc and MW models, it is likely that other processes 
(shock-heating by SNe and stellar winds, photo-ionization, mag- 
netic fields) can play a significant role in shaping the gas disper- 
sions and density distribution. We therefore defer a more detailed 
comparison with these observations until the models include some 
of these effects. However, preliminary experiments show that while 
the systematic values discussed above can shift, the qualitative con- 
clusions remain intact as other feedback mechanisms are intro- 
duced. 

4 DEPENDENCE ON MODEL PARAMETERS 

In this section we show that the results summarized in ^do not 
depend sensitively on the assumed local star formation law ( ^4.1^ 
or the precise feedback parameters adopted ( ^4.2^ . We focus on 
the star formation history when presenting these results. We again 
focus on the HiZ model since it is the most self-gravitating and 
therefore its SFH tends to be the most sensitive to variations in the 
simulation parameters. However, we carried the same experiments 



" In detail, we find there is a nanow core in the distribution and a broader- 
than-lognormal wing; the distribution is better tit by an exponential at high 
columns, with PilogNn) oc exp [— | logA'^/A'o 1/0.22]. This is broadly con- 
sistent with randomly distributed "patchy" obscuration within clouds. 



for the MW-like simulation and found comparable results, which 
are also shown below. In the appendix we show that our results 
also do not depend strongly on how we numerically implement the 
stellar feedback. 



4.1 Dependence on the Local SF Law 



Figures [6|7[ shows how the star formation history in feedback- 
regulated simulations depends on the local star formation prescrip- 
tion used at high densities. For our fiducial rjp = r],, = \ model. 
Figures [6|7[ vary the star formation efficiency in dense gas e, the 
power-law slope of the star formation law, and the threshold den- 
sity for star formation no (eq.[2](. 

The key result in Figures[6j7|is that there is very little depen- 
dence of the SFH on the high density star formation law. Specifi- 
cally, Figure [6|(Ze/!' panel) shows results for our canonical value of 
e = 1.5%, a larger value of 6%, and a smaller value of e = 0.35% 
(we have also examined several intermediate values). This range of 
e corresponds to a factor of 20 different star formation timescale at 
fixed density. We find, however, at most ~ 30% differences in the 
SFR once the system has reached approximate equilibrium. In the 
MW-like model (Figure |7| the conclusion is identical. Secondly, 
we vary the power-law index of the local SF law {middle panel). 
In our canonical implementation, p» oc p/fdyn oc p' '; we compare 
this to simulations with p, oc p/to oc p'" and p, oc p^", normalized 
such that p, is identical at the threshold density po. There are early- 
time differences in the star formation histories, but given the magni- 
tude of the change to the star formation prescription the results are 
broadly similar. The biggest change appears when p, oc p, i.e. when 
the gas consumption timescale is constant, independent of density; 
in this regime, the gas cannot necessarily be consumed quickly on 
small scales, so the collapse from large to small scales is no longer 
the dominant rate-limiting step in star formation (a slightly larger 
exponent, e.g. p, oc p' much more closely resembles the canoni- 
cal oc p' case). We show this for the MW-like model in Figure|7] 
comparing p* oc p''^, p* oc p^", p* oc p' The relative differences 
in the p, oc p' " are even smaller, and making the exponent just 
slightly super-linear (oc p' ') gives a nearly identical SFR to oc p'^. 
Finally, we vary the SF density threshold no {right panel); from 
our canonical value of lOOcm^'', we also consider no — 25cm~^ 
and no = 2500 cm^^^ (with other intermediate values sampled as 
well). At early times, before the initial conditions have been re- 
placed by a self-consistent equilibrium, the SFR is higher with a 
lower threshold (unsurprisingly). However, once enough time has 
elapsed for gas to collapse to high densities and initiate significant 
feedback, the SFHs are again nearly identical despite a factor of 16 
change in the threshold for star formation. The same result obtains 
in the MW model, for which (given the lower mean density) we 
vary no = 10 — 1000. Moreover, for a MW-like model with similar 
resolution, [Saitoh et al.[p008l > find the same result in a more limited 
study varying the small-scale star formation efficiency, but with a 
different simulation algorithm and different feedback mechanism 
(SNe) implemented. We have also repeated these experiments for 
different values of r]p (= 1/3, 4, 10) and r]y {= 2 and continuous 
acceleration), and reach similar conclusions in each case. 

Our interpretation is that the weak dependence of the global 
star formation rate on the small scale star formation model is a 
consequence of the turbulence driven by stellar feedback, and the 
self-regulation to 2 ~ 1 (see, e.g., [Thompson et al.|2005| . Specif- 
ically, gravity causes gas to collapse to high density, where some 
of it forms stars, while most of the gas is driven back out to lower 
densities by feedback. The key step that regulates the star formation 
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Figure 8. Cumulative gas mass fraction above a given density n for the (liigh-resolution) HiZ model for variations in the small-scale star formation law (Fig. [6} 
and feedback efficiency (Fig.|9|. Left: Density distribution for different values of the star formation efficiency e: for smaller (larger) e, more (less) mass must 
collapse to high densities for the star formation to self-regulate (the high-p cutoff is set by resolution limits). Middle: Density distribution for different values 
of the momentum deposition per unit star formation (?;,,; eq.|5j. For larger rjp, gas is more efficiently removed from dense regions. Right: Density distribution 
for different values of the threshold density for star formation no ■ Larger hq requires that the gas collapse to somewhat higher densities before the star formation 
can self-regulate. 



rate is this cycle of collapse and expulsion, which has a timescale 
~ the global dynamical time of the galaxy - this is also the decay 
timescale for large-scale turbulence in the galaxy. The details of 
feedback on small scales should also not be important, so long as 
it is sufficient to self-regulate (compare our result to [Saitoh et al.| 
|2008[ >. So long as the star formation timescale at the threshold den- 
sity is small compared to the global dynamical time (i.e., e not too 
small and po > the (p) of the galaxy) and the threshold is well- 
resolved numerically (i.e., po is not too large), the star formation 
rate is insensitive to the details of the small-scale star formation 
law. 

More generally, if the support needed to maintain stability 
against runaway star formation is set by the luminosity/mass in 
young stars, the SFR can self-regulate to 2 ~ 1 . For example, if the 
SFR set by the small-scale physics is too low to maintain Q ^ \ 
given the large-scale conditions, gas simply collapses further to 
slightly higher densities until the required feedback power is gener- 
ated, sufficient to halt further collapse. The high density star forma- 
tion law thus determines some of the properties of the high density 
gas, but not the global SFR. 

Figure [8] supports this interpretation by showing the cumula- 
tive gas density distribution (mass fraction > n) for different values 
of £ (left panel), no (right panel), and the feedback parameter rjp dis- 
cussed in the next section (middle panel). Figure[8]shows that when 
the high-density star formation efficiency e is smaller (larger), the 
gas distribution adjusts so that there is more (less) mass at high 
densities, so as to produce a similar total SFR (as in Fig.|6]) When 
the threshold density no is varied, the mass at high densities shifts 
accordingly. For example, increasing no causes the gas that would 
have formed stars at the previous threshold to collapse to somewhat 
higher densities before it begins to form stars. 

Note that the mass at low densities is nearly unchanged - the 
disks are not in global collapse (they are regulated by feedback), 
but the gas locally collapses to the densities needed to maintain the 
same SFR. For this reason, the Schmidt law predicted by each of 
the models in Figure[6]is nearly identical. They have the same range 
in surface densities (set by the initial conditions and exhaustion via 
star formation, which must be the same since they have the same 
star formation history), and so self-regulate at the same SFR. 

[Schaye et al.| ( [20T0l l, using much lower resolution cosmolog- 
ical simulations, also find a galaxy wide star formation rate that 
is independent of the details of the small scale star formation law 



employed. However, in their case, because star formation laws are 
applied globally (on >kpc scales), it is the global gas mass that 
self-adjusts (e.g. lowering the star formation efficiency leads to in- 
flows larger than the SFR building up the global gas mass until the 
SFR is similar to the cosmological inflow rate), so the systems do 
not necessarily obey the observed Schmidt-Kennicutt relation. In 
our case, neither the SFR nor global gas mass varies; what does al- 
ter is the gas fraction at the very highest densities available to the 
simulations. 

4.2 Dependence on the Feedback Efficiency 

Figure |9] (Figure [T0| shows how the star formation history of our 
HiZ (MW) model depends on the feedback parameters ry,, and 
rj,. (eqs. |5] & |7j and on whether we implement the momentum- 
feedback continuously or via kicks (left and middle panels). All of 
the variations are with respect to our standard rjp — rjv — I model. 

Figure |9] (left panel) shows that simulations with rjy = 1 and 
rjy = 2 produce very similar SFHs. These two simulations both have 
rjp — I and thus have the same momentum injection rate. Physi- 
cally, the similarity in their SFHs arises because the particles in- 
teract with ambient gas and share their momentum efficiently. The 
end result is that clumps being destroyed by stellar feedback have 
velocities comparable to the escape velocity from the clump, rela- 
tively independent of the initial velocities we input. Figure |9] also 
shows a comparison of two different methods of implementing the 
same momentum flux: particle 'kicks' and continuous acceleration 
(see ^^23). It is reassuring that these two methods produce quite 
similar results - this again highlights that the critical parameter is 
the rate at which momentum is deposited into the ISM, not pre- 
cisely how it is deposited. The MW-like model in Figure [TO] gives 
identical conclusions (the dependence is even weaker in this case). 

Figure|9](m!iirf/e panel) also compares simulations with varied 
momentum-injection per unit star formation rjp, from rjp — 0.3 — 10 
(all at fixed rj^ — 1). As expected, the quasi-steady SFR decreases 
as the efficiency of momentum-injection increases. However, the 
decrease in the SFR is rather mild, with a factor of ~ 2 change 
in SFR over a factor of 10 in rjp. The same scaling holds in the 
MW-like model in Figure fTO| N aively one m ight expect an inverse 
linear scaling M, oc rjp' ( ^Thompson et al.|2005| ). Specifically, the 
turbulent energy dissipation rate per unit area is 

-^~{E)5v'n~G{E)'5v (11) 
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Figure 9. Star formation rate as a function of time for the HiZ model for different feedback parameters and resolution. Left: Variations in the initial kick 
velocities at fixed momentum-loading. We also show a simulation with stochastic particle "kicks" replaced by continuous acceleration of all particles. In all 
of these simulations, the gas shares its momentum with the rest of the surrounding clump and thus produces similar dynamics. Center: Variations in the 
momentum deposition per unit star formation (r;,,; eq.|5j. The star formation rate decreases by less than a factor of 2 over a factor of 10 in rjp. Right: Variations 
in resolution. Increasing the particle number from ~ 10^ to ~ 10^ (our "intermediate" vs "high" resolution) increases the star formation rate at early times by 
a moderate amount (~ 20 — 40%). But after about one orbital time, the results are quite similar. Going to yet higher resolution gives nearly identical results. 
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Figure 10. As Figure[9] but for our MW model. Dependence on rj,, and res- 
olution is even smaller than the HiZ case; the dependence on rjp is similar. 



and the optical depth in those regions - to decrease roughly as ri~ . 
This is why both the momentum input cx rjpTiR and the star forma- 
tion rate oc (»7;,riR)~' have only a weak dependence on -qp. This 
property of our numerical simulations is one of the most significant 
differences between our results and previous analytic treatments 
of star formation regulated by radiation pressure ( [Thompson et al.| 
|2005[ l. It is important to stress that this self-regulation to achieve 
the same star formation rate relatively independent of the feedback 
parameter r/p is only a property of models in which the momen- 
tum injection rate is proportional to the gas surface density (eq.|5]l; 
that is, it is only a property of feedback by radiation pressure, not 
momentum injection associated with supemovae or stellar winds. 



where (E) is the mean surface density of the disk and we have 
assumed 2 ~ 1 and that turbulence dissipates on a crossing time 
~ /i/(5v ~ r2~'. The total momentum injection rate scales with the 
SFR as P ~ (1 +ripTiR)L/c ~ (1 + r/p Ekir) e,M* c, where e, ^ 
4 X 10""* — L/Mt, c^. If the momentum-injection ultimately drives 
turbulent motions with random velocity 5v, the associated energy 
injection rate is ~ PSv. Balancing this against the dissipation in 
equation [TT] we find that the star formation rate per unit area Esfr 
is given by 



4.3 Dependence on Resolution 



G{E)^ 



,c{l+ripTiB,) ' 



(12) 



For parameters relevant to Figure [9] »7priR > 1 and so equation 



implies that M, oc (riR?7p)~'. This does not, however, imply 
M* oc 77,7'. One reason is that equation 1 1 2| neglects the possibil- 
ity of significant cancellation in colliding/canceling flows. More 
importantly, however, for (say) larger rjp, the fraction of mass at 
high densities decreases because feedback is more effective. This 
is shown explicitly in Figure [8] (mirfd/e panel): as rip increases, the 
density distribution cuts off more sharply at high n. As a result, the 
optical depth tir in the regions of massive star formation decreases. 
This demonstrates that the momentum input oc rjpTm, and thus the 
star formation rate, must scale sub-linearly with rjp (as in Fig.|9]l. 
Assuming that feedback removes gas from high to low density at a 
rate oc rjp, we. would expect the fraction of mass at high densities - 



Figures [9| 10| (right panel) shows how the star formation histories 
of our fiducial rjp = rjy— I HiZ and MW models depend on particle 
number, with Np^ = 2 x 10*, 2 x lO', and 2 x 10**. The basic evolu- 
tion of the SFR is very similar in all cases. The SFR is ~ 25 — 40% 
higher at early times in the N^m = 2 x lO' simulation relative to 
Apart = 2 X 10*, but there is a much smaller change going to yet 
higher resolution. Moreover, after a few dynamical times, all of the 
simulations have a comparable SFR. We find similar convergence 
for different galaxy models and different feedback parameters. 



5 THE GLOBAL SCHMIDT-KENNICUTT LAW 

Figure [TT] compares the global Kennicutt-Schmidt law predicted 
by simulations with and without our stellar feedback model. We 
measure Esfr = A/, (< Rsfr)/'KRlf,- as a function of Egas = Mgas(< 
/?sfr)/7!"^sfi, where /?sfr at each time is defined as the half-SFR ra- 
dius viaM, (< Rsh ) =M,/2. This radius is chosen to loosely corre- 
spond to the half-optical or half-Ha radii used in various observa- 
tional studies, but adopting a different choice (e.g. the half-gas mass 
radius) primarily shifts the models along the relation. The numeri- 
cal results are shown every Myr. The numerical results in FigurefTT] 
are compared with several different observational inferences of the 
K-S relation: the best-fit power-law relations from low-redshift in 
|Kennicutt| ( (l998| l and high-redshift in |Bouche et al.| ( |2007[ >; |Genzel| 
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Figure 11. The global Kennicutt-Schmidt relation between star formation rate density and gas surface density in our simulations. Left: Without feedback. 
Right: Simulations with our feedback model for a range of parameters (see Table |2j. Each row is a different galaxy model (see Table [TJ. In each panel, a 
point corresponds to one time snapshot in the simulation, evenly spaced in ^ 10* yr intervals (starting after two dynamical times). The surface densities are as 
viewed face-on, averaged within the circular radius that encloses 1/2 of the star formation. Solid lines show the fits to the data in |KennicuttH 1 998| and updated 
with high-redshift galaxies by |Genzel et al.|j2010) . Grey shaded region shows the 90% completeness range at each gas surface density from the compilation 
of the systems observed in those two works as well as the compilations in |Bigiel et al.|j2008) and |Daddi et al.|j2010} . Without feedback, the gas experiences 
runaway collapse and is consumed in less than a dynamical time, predicting star formation rate densities in excess of the observed KS-relation by factors of 
~ 100 — 10*. With feedback, the gas disks quickly self-regulate and reach an approximate equilibrium comparable to that observed. 



|et al.|pOTol l, together with the 10 — 90% interval of all points from 
the combined compilations in those studies as well as |Bigiel et al.| 

'* Note that the shaded range falls below the best-fit power-laws at low 
surface densities because the power-law fits did not study the low surface- 
density 'cutoff due to a low molecular fraction. 
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with and without feedback have identical initial conditions in each 
case. 

Without feedback our simulations predict a SFR surface den- 
sity well in excess of the K-S law (see also Fig.[4]l. Absent feed- 
back, the gas cools and cannot avoid runaway collapse; most of the 
gas is consumed into stars in a single dynamical time, leading to 
SFR surface densities ~ 10 — 10* times larger than observed. 

In contrast, the simulations with feedback lie close to the ob- 
served relation at essentially all times. This is true over the range 
of feedback parameters and resolution we have studied; the runs in 
Figure[TT]span a range in JTp = 0.5 — 10, 77,, = 1 — 4, star formation 
law variations as in Figure Im and resolution (A'paiticies ~ lO' — 10^). 
Varying the simulation parameters for each galaxy model tends to 
shift the systems along the KS relation, rather than dramatically off 
the relation. For each galaxy model, there is a broad dynamic range 
in S covered; in particular, the high-z simulations lie on the ob- 
served relation over multiple decades in density. The average slope 
of the relation (if we consider all galaxies together) is quite similar 
to that observed; however we see that there can be significant vari- 
ation in th e slope within galaxies, also commonly observed (see 
Bigiel et al. 20081.'^ Altogether, the feedback-regulated simula- 
tions lie on the observed KS relation over a dynamic range from 
~ Egas ~ 10' - IO'^Mq kpc"-. The scatter about the Schmidt law 
predicted is also similar to that observed, about 0.5 dex. 

6 DISCUSSION 

We have presented a new numerical method for treating stellar 
feedback in hydrodynamic simulations of galaxies. We have imple- 
mented this method in the SPH code Gadget-3 but our approach is 
general and can be utilized in both Lagrangian and Eulerian codes 
(see Our stellar feedback model is motivated by the physics of 
feedback in dense environments: under these conditions, gas cools 
rapidly and the primary dynamical influence of stellar winds, su- 
pemovae, and the stellar radiation field is the momentum they im- 
part to the ISM. In addition to formulating the general method, we 
have carried out a detailed study of the properties of this stellar 
feedback model in isolated (non-cosmological) disk galaxy simu- 
lations, from models motivated by massive z ~ 2 galaxies forming 
stars at ~ 100 — 300 Mq yr~ ' to models of SMC-like dwarf galax- 
ies. These disk galaxy calculations are not intended to be quantita- 
tively applicable to real systems; rather, they illustrate our method 
and demonstrate the critical importance of including stellar feed- 
back by momentum injection. In a future paper, we will combine 
the method developed in this paper with more standard treatments 
of supernova and stellar wind heating, to produce a more compre- 
hensive stellar feedback model. 

High-resolution numerical simulations of isolated galaxies 
and galaxy mergers, as well as cosmological "zoom-in" simula- 
tions, can readily resolve the formation of numerous dense gaseous 
clumps via gravitational instability (provided the cooling to low 
temperatures ~ 100 K is not artificially suppressed); we dub these 
clumps GMCs although we do not include the physics of molecule 
formation in our simulations. In observations of nearby galaxies, 
most of the star formation occurs in GMCs - this is also true in our 



For example, the MW-Iike simulation is significantly steeper at low den- 
sities than the median relation, but quite similar to spatially resolved obser- 
vations in M51 (Kenn icutt et al.|2007^ . However, this is the regime where 
we expect other physics (e.g. the addition of SNe feedback and possibly the 
effects of detailed molecular chemistry) may become important. 



simulations - and thus it is important to have at least an approxi- 
mate model for GMC disruption by stellar feedback. 

To model stellar feedback, we implement an on-the-fly cluinp 
finding algorithm to identify high density star-forming clumps (i.e., 
GMCs). We then deposit momentum into the surrounding gas at a 
rate proportional to the radiation produced by young stars in the 
clump; this force is directed radially away from the center of the 
GMC. More precisely, the force we apply scales as ~ tirL/c (eq.jsjl 
where tir is the optical depth of the clump to IR photons and the 
stellar luminosity L is calculated as a function of time given the 
stellar ages using Starburst99 models. Although our model is quan- 
titatively motivated by radiation pressure on dust, the momentum 
flux from supernovae and massive stellar winds can also be sig- 
nificant. We will study the relative importance of these different 
feedback mechanisms in detail in a future paper. 

The model we have developed is distinct from the stellar feed- 
back models used in most of the galaxy formation literature. First, 
we input momentum, rather than thermal energy, into the ISM 
around young stars. The motivation for this choice is that mo- 
mentum, not energy, is the relevant conserved quantity in dense, 
rapidly-cooling gas. Moreover, the feedback we implement scales 
with the local surface density of the GMCs, as expected for the ra- 
diation pressure produced by stellar photons as they are degraded 
by dust from the UV to the far IR (eq.|5j. As summarized below 
and in ^ we find that this surface density dependence is critical to 
the evolution of our galaxy models. 

In our study of stellar feedback in isolated galaxies, we do 
not "turn off" hydrodynamic forces, cooling, star formation, and/or 
other physics in the gas to which the feedback is applied. By con- 
trast, many stellar feedback implementations in the galaxy forma- 
tion literature either turn off hydrodynamic forces in winds for 
some free-streaming length (typically such that winds escape the 
galaxy) or turn off cooling and star formation in supernova-heated 
gas for some period of time. In such models, the induced veloci- 
ties on galactic scales are essentially determined by hand (through 
adjusting the relevant parameters) as is the presence/absence of a 
global galactic wind driven by stellar feedback. In our model, the 
single key parameter is the momentum supplied to high density 
gas around star clusters ~ the resulting galaxy-wide turbulence, the 
properties of galactic winds, etc. are all predictions of our model. 
The model remains "sub-grid," but on the scale of individual molec- 
ular clouds rather than the galaxy as a whole. 

We are able to directly model the stellar feedback without arti- 
ficially modifying the underlying equations for several reasons. The 
high resolution in our simulations allows us to partially resolve the 
multi-phase ISM structure: since star formation is spatially inho- 
mogeneous, the stellar feedback is as well, which self-consistently 
maintains a turbulent and multi-phase ISM structure (Fig.[TJ. Per- 
haps more importantly, the feedback is momentum-driven and the 
forces are directed away from the centers of local gas overdensities 
(GMCs), the sites of massive star formation. As a result, the feed- 
back is effective even in dense regions of the ISM, in which the 
cooling time is much shorter than the dynamical time. In standard 
treatments of feedback by supernovae, the feedback is inefficient 
in dense regions because the thermal energy supplied by super- 
novae is rapidly r adiated away ([Thack er & Couchman'2000' 'Gov-| 
|emato e t al.||2007l|Ceverino & Klypin 2009; Boumaud et al. 2010} 
Teyssier et al. 2010: Brook et al. 2011). Thus many simulations that 
nominally include stellar feedback do not in fact have feedback that 
is quantitatively of the correct order of magnitude. 

As a first assessment of the implications of our new stellar 
feedback model, we have used it to study star formation in a wide 
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range of (non-cosmological) disk galaxy models, representing sys- 
tems ranging from SMC-like dwarfs, to the MW and local LIRGs, 
through to massive high-redshift gas-rich disks. The disks are ini- 
tially pressure-supported, but cool rapidly to < lOOK and collapse 
into a wide spectrum of GMCs. Absent stellar feedback, we find 
that GMCs undergo a runaway gravitational collapse to high den- 
sity; star formation proceeds on approximately a single galaxy av- 
eraged dynamical time (Fig. |4]l, a result that is dramatically in- 
consistent with observations (Fig.[TTJ. However, with feedback in- 
cluded, the GMCs dissociate once a modest fraction of their mass 
has turned into stars and the galaxy develops a turbulent, multi- 
phase, ISM as long as gas remains. Quantitatively, the turbulence 
in the ISM maintains marginal stability to self-gravity, i.e., 2 ~ 1 
(Fig.|5j. Moreover, the galaxies self-regulate and approach a quasi- 
steady state star formation rate that is consistent with the observed 
Kennicutt-Schmidt relation over a dynamic range of several orders 
of magnitude in surface density (Fig.|l l|l. 

Our numerical results are reasonably consistent with the ob- 
served global Kennicutt-Schmidt relation nearly independent of the 
high-density star formation law used in the simulation (Fig. [6] & 

This is in contrast to many results in the literature, where free 
parameters in the high-density star formation law are adjusted to 
approximately reproduce the Kennicutt-Schmi dt relation ([Springel 



& Hemquist 2003a Govemato et al.|[2004[ [Dubois & Teyssier 
2008 1 |Agertz et al.|2009| . A weak dependence of the star forma 



tion rate on the high-density star formation law is important for 
developing a predictive galaxy formation model. It is otherwise dif- 
ficult to disentangle results that are due to the physics of star forma- 
tion and/or feedback from those that are due to particular numerical 
choices/parameters. 

Our star formation model is that gas turns into stars in dense 
regions above some threshold density po at a rate p, oc p" . Varying 
the normalization of this relation (the high-density star formation 
efficiency) by a factor of ~ 20, varying po by a factor of ~ 100, 
and varying the power-law index n in the range 1—2 changes the 
quasi-steady state star formation rate by < 50% (Fig. [6]l (this of 
course requires that the threshold density be well-resolved). Physi- 
cally, this weak dependence arises because the condition for quasi- 
steady state star formation is that the momentum injection rate by 
stellar feedback is sufficient to maintain the ISM at 2 ~ 1 . Reach- 
ing 2 ~ 1 requires a particular turbulent velocity 5v, and thus a par- 
ticular momentum injection rate, for a given set of global galaxy 
properties (eq. [Toj. Variations in the high-density star formation 
law are compensated for by slightly more or less gas collapsing to 
high densities (and differences in how dense the gas becomes be- 
fore GMCs are dissociated), so as to produce the same momentum 
injection rate and hence the same global star formation rate (Fig.jSj. 

The key parameter that determines the efficacy of the stel- 
lar feedback in our model is the normalization of the momentum 
injection rate, 77,, (eq. [sjl, where p ~ (1 + ?7p'nR)L/c; physically, 
77p < 1 corresponds to photons leaking out of regions with lower- 
than-average surface densities while 77,, > 1 corresponds to the ef- 
fects of additional momentum sources (e.g., supernovae and stellar 
winds) and/or insufficient resolution of the highest optical depth 
(tir) regions. Numerically, we find less than a factor of 2 change 
in star formation rate over a factor of ~ 10 in rip (Fig. |9j. Physi- 
cally, this is again because maintaining 2 ~ 1 requires a particular 
momentum injection rate and thus a particular star formation rate. 
Variations in are compensated for by the surface densities and 
thus optical depths tir reached in dense star-forming clumps, main- 
taining approximately the same momentum injection rate indepen- 
dent of rip (Fig.|4]&[|]l. 



Although we have emphasized the importance of momentum 
input throughout this paper, this is clearly only part of the impact 
of massive star formation on the ISM of galaxies. Which feedback 
process is the most important depends on the galaxy mass, gas frac- 
tion, etc., and on the specific science question of interest. Heating 
by photoionization and supernovae, and their effect on molecule 
formation, are critical physics to include in the formation of the first 
stars as well in studies of lower-density gas characteristic of dwarf 
galaxies and the outer parts of more massive disks. In the diffuse 
ISM and the halos of massive galaxies, additional pressure sup- 
port from cosmic rays and/or magnetic fields may also be impor- 
tant. The model presented in this paper is most directly applicable 
to dense gas in the central kpcs of massive, enriched and evolved 
systems, in which cooling times are short and molecular fractions 
are order unity. And even in these regions, the model here under- 
predicts the temperatures of the "hot" diffuse ISM {T > 10^ K); 
this gas is likely to be heated by shocks from SNe explosions and 
fast stellar winds, with v ~ lOOOkms^'. Although at a given in- 
stant, this phase represents only ~ 1% of the gas mass, it can have 
important effects on the generation of galactic super-winds. In a 
subsequent paper we will study the combined effect of stellar ra- 
diation, stellar winds, and supernovae, with the goal of developing 
a more widely applicable stellar feedback model for use in galaxy 
formation. To extend the study here from idealized disks to disks 
forming over cosmological timescales, it will also be important to 
incorporate the the realistic cosmological effects of gaseous halos 
and cold-flow accretion as well as galaxy mergers. 
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APPENDIX A: ADDITIONAL NUMERICAL TESTS 

In this appendix, we discuss a number of additional numerical tests 
performed to ensure that our conclusions are not sensitive to the 
precise implementation of the feedback algorithm. 
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Figure Al. Example star formation histories for the HiZ galaxy model for 
variations in the details of our numerical method. 7b/?; Otherwise identical 
simulations with different choices of N^mi (linking length for clump finder), 
A'* (number of lengths to smooth/cut the radiation field), and Wjeans (numer- 
ical pressure support to prevent artificial fragmentation). Bottom: Otherwise 
identical simulations with variations in the direction in which the feedback 
is applied. Applying the force radially outwards from the nearest gas den- 
sity peak, stellar center of mass, or center of light give nearly identical re- 
sults. The "isotropic" models choose the direction of each kick randomly. 
When individual kicks are relatively rare and large (77,- > 1), this is some- 
what less efficient than radial kicks but still stirs turbulence and slows down 
star formation. By contrast, in the limit of continuous acceleration, random 
isotropic forces cancel, impart little net momentum, and the star formation 
rate is only a factor of ~ 2 smaller than that found in the absence of feed- 
back (Fig.|4). 

Figure [ATj shows examples in which we modified the friends- 
of-friends search used to identify the nearby density peak that is 
the origin for our feedback (see §2.2. Specifically, we varied the 
parameter A'smi between 1 — 5 (our standard model adopts Wsmi = 3) 
- this defines the number of smoothing lengths in which to search 
for a more-dense gas particle in each iteration. Within this range we 
do not see significant differences in most of the identified density 
peaks; nor is there a significant effect on the star formation history 
(Fig. [ATJ. We have also modified the resolution-dependent pres- 
sure floor used to avoid artificial collapse, through the parameter 
Means which represents the minimum number of smoothing lengths 
which resolve the Jeans length (eq.[T]l. Because all of our feedback- 
regulated runs have large feedback-induced random velocities, this 
pressure floor makes no significant difference in these cases (see 
FigjAlJ. It does, however, determine the smallest scale of resolved 
structure in simulations without feedback (see e.g. [Robertson &| 
|Kravtsov|2008l >. 



To reduce the code run-time and ensure that feedback is rela- 
tively local to massive stars, we only apply the feedback to gas par- 
ticles within Nt =3 smoothing lengths of any star particle ( §2.2.2| l. 
We have experimented with this value in the range 2 — 20 (the latter 
value including almost all of the gas). Because the optical depths 
and momentum deposition are dominated by the gas closest to the 
stars, this choice has very small effects on the ISM properties and 
star formation history (Fig. |Al[ ). 

One aspect of our method that can have a more significant 
influence on the results is the direction in which particles are ac- 
celerated when the feedback is applied. In the standard model, we 
accelerate particles away from the gas density peak identified in 
the friends-of-friends search ( §2.2.2^ . We have experimented with 
changing the origin for this force to be the center of mass of the 
stars or gas or the center of luminosity of the star-forming clump. 
We have also considered models in which the force vector is ori- 
ented along the local density gradient (appropriate in principle for 
arbitrary geometries, in the regime where tir falls below unity 
along that gradient). We find that these models give nearly indis- 
tinguishable results, as is shown in Figure [AT| This is essentially 
because most of the stars are concentrated near the clump center 
(by any of these metrics) for most of the time when feedback is 
important. On small scales in galaxies (or interior to GMCs), how- 
ever, the dynamical time can be much shorter than the lifetime of 
massive stars, and so it is possible that large separations could arise 
between massive stars and the gas from which they formed. In this 
case it would be better to determine the direction of the force using 
the local peak in stellar luminosity. 

We have also considered experiments in which we com- 
pletely ignore the clump density information and kick particles with 
isotropic, random directions (rather than away from clump centers). 
In our standard model (77,. = 1), kicks are somewhat rare but have 
large initial velocities, so the coherent momentum imparted with 
each kick is still large (even if it is randomly directed). In this case, 
the feedback is somewhat less effective than in our standard model 
and so the star formation rate is somewhat larger (Fig. |ATJ. If the 
individual kicks particles receive are much smaller (but more fre- 
quent), the coherent momentum imparted will be reduced if each is 
independently randomly oriented. As a result, in the limit in which 
we continuously accelerate particles in random directions (rather 
than imparting discrete kicks), we find that the feedback has little 
effect. The star formation history is similar to models with no feed- 
back (Fig. |AI^ . This highlights the importance of properly applying 
the feedback radially away from the center of mass/luminosity of 
massive star clusters. 

There are other aspects of our simulations that are uncertain, 
independent of the stellar feedback model. For example, our cool- 
ing function at low temperatures is not exact, since we do not ex- 
plicitly follow chemical networks. We have therefore considered 
various arbitrary changes to the cooling function: setting A(r) be- 
low 10"* K to a constant median value, or simply forcing all gas at 
high densities to a minimum temperature ~ 100 K. These intro- 
duce < 20% changes in e.g. the star formation history, since in all 
cases the cooling time is short relative to the dynamical time. On 
the other hand, removing fine-structure cooling entirely (effectively 
producing a cooling floor at 10"* K) dramatically changes the behav- 
ior in the MW, Sbc, and SMC cases, since this temperature floor is 
sufficient to artificially prevent collapse to high densities. However 
it makes little difference in the HiZ case because the requirement 
fori2> lisc,> 30-50kms"'. 
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Figure Bl. Modified coupling of photon momentum in a self-similar, 
"leaky" medium with a broad density distribution. The true "boost" to the 
coupling Teff defined such that P = T^ifL/c, is plotted as a function of the 
mean tq, for a source at the center of a medium with a random distribution 
of densities that obeys a power-law (PL; Eq. |Bl) or lognormal (LN; Eq. |B3^ 
PDF with logarithmic dispersion a. Dotted line shows r^ff = tq, the expec- 
tation for a completely homogenous medium. For cr < 1 in the lognormal 
model, or < 0.5 in the power-law model, r^ff <x tq (with relatively small 
normalization corrections comparable to small variations in our rjp parame- 
ter). At larger dispersion, the scaling becomes sub-linear, with r^ff oc t^^" 

(PL) or Teff oc .^^"^o/4ct ^lN). The dispersion in ultra high-res simulations 
(and observations) corresponds to cr 0.5. 
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Figure B2. Effects of IR photon "leakage" as calculated in Figure [BT] on the 
SFH of simulated galaxies (here the HiZ model). We compare our standard 
model (Teff = ro) to models using the calculated Teff (ro) from Fig. |Bl| for 
an assumed universal <j and functional form of the density distribution as 
in that Figure. We also compare a model with the power-law/exponential 
distribution (Eqn. |Bl| and a calculated on-the-fly as the dispersion in ln(p) 
within the identified clump radii Rclump. Larger dispersion in columns leads 
to more leakage and higher SFRs, but the effect is weak and, in the regime 
of interest here, similar to a choice of rjp slightly somewhat less than unity. 



APPENDIX B: THE EFFECTS OF PHOTON "LEAKAGE" 

A subtle complication in modeling the effects of radiation pres- 
sure arises if the ISM is truly inhomogeneous on all scales, includ- 
ing those well below what we model. A patch of ISM that appears 
smooth in the simulations, with some average optical depth to, may 
more accurately (at infinite resolution) exhibit a distribution of lo- 
cal columns, including some optically thin lines of sight that could. 



in principle, allow photons to "leak out" at a rate much higher than 
the nominal exp(— to) expectation. This would potentially lower 
the effective "boost" to the radiation pressure from tqL/c to some 
TeffL/c. It is straightforward to show that leakage will not signifi- 
cantly change the total energy absorbed and therefore the IR lumi- 
nosity density - once ro is large, it is generically true for essentially 
any realistic distribution of optical depths that most of the incident 
opticalAJV radiation is initially absorbed (whether the escape frac- 
tion is a few tens of percent or vanishingly small makes only a tens 
of percent change to Lir). The concern is rather that IR photons 
will tend to escape along optically thin sightlines before they scat- 
ter sufficiently numerous times to impart the full ~ to momentum 
boost. |KrumhoIz & Matzner| ( [2009| ) for example argue that the lat- 
ter effect means that the effective Teff can never be larger than a 
few, even when to ^ 1. But they do so by assuming an order-unity 
fraction of un-obscured sightlines, independent of to. We therefore 
consider this effect in more detail in this Appendix. 

As discussed in the text, the case of a perfectly homogenous 
density distribution with a source at the center is trivial. The opacity 
along all sightlines is to, so the photon scatters an average ofN — Tq 
times as it performs a random walk to diffuse out of the sphere. For 
a random walk, the net momentum flux directed radially away from 
the central source is just TeffL/c where Teff — Vn, which in this case 

gives Teff = To. 

But now consider a case with an inhomogenous density dis- 
tribution. Unfortunately in general, calculating the precise Teff for 
any inhomogenous density distribution is complex and cannot be 
solved analytically - it requires a full radiative transfer solution 
for each specific density distribution. However, we can make con- 
siderable progress and obtain reasonably general scaling laws with 
some simplifying assumptions. Consider a "cloud" of ISM with a 
well-defined mean to (which we can measure easily in the sim- 
ulations) enclosing a source at its center; for convenience (with 
no loss of generality) define the cloud radius to be ^o = 1- Now 
define the "true" distribution of optical depths across all sight- 
lines within the cloud to be dP{T\To). Finally, assume that the 
cloud is self-similar with structure on all scales. In this limit, the 
distribution of local densities or At/AI (where I is the distance 
a long a line of sight) is just dP{dT/Atj = ^/'(t/^o j to/^o) = 
dP{T\To). We can perform the following relatively simple calcu- 
lation: take a population of photons starting at the center, with ini- 
tial random directions. For each, draw a random At /Al (equiva- 
lently, line-of-sight-averaged density), and determine in standard 
monte carlo fashion the distance the photon travels before being 
absorbed (for a uniform random variable p between zero and unity, 
Al = — In (1 — p)/[dT/d^]). At each "scattering," determine a new 
random direction for the photon to be re-emitted, and record the lo- 
cally coupled momentum as the negative of the change in the pho- 
ton momentum. This is iterated until all photons escape the sphere; 
using a large monte-carlo sample of "photons," then, we can deter- 
mine quantities such as the average number of scatterings (W) and 
the average momentum imparted (or effective boost Teff).'* 

Of course, we still need to specify some distribution 
dP{T I To). Fortunately, we can make a reasonable estimate: in ultra- 
high resolution simulations, we can calculate for example the form 
of dP{T I To) for each molecular cloud in the simulation, using a 
large number (~ 1000) lines-of-sight and integrating the simulation 



" Specifically we are interested in the net momentum imparted radially 
away from the center. As expected, all other components of the coupled 
momentum average to zero. 
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column along each sightline. In § |2.2.2| we discuss this process and 
note that the resulting column density distribution on a per-cloud 
basis can be well-approximated by a near-universal function 

./'(.|.c).J-exp(-MlM)^^I (Bl) 
2(T V a / T 

with a (the standard deviation) ranging from 0.25 — 1.0 (0.1 — 
0.4 dex) with a median a — 0.5 (0.22 dex). This is very similar 
to the distribution of columns estimated in much higher resolution 
simulations of individual GMCs and sub-cloud clumps (often with 
very different physics included; see e.g. |Ostriker et al.|200T) , and 
to observational estimates of the column density distribution across 
observed GMCs ^Wong et al.|2008i [Goodman et al.|2009l l. 

Note that this distribution is exponential in log(T) - really, 
the key behavior is that the number sightlines at small or large r 
falls off as a power-law in (t/tq), rather than an exponential or 
log-normal. This is important, as we will see below. 

Given some a, then, it is straightforward to perform the Monte 
Carlo calculation of Tcff described above. Figure [BT| shows the cou- 
pled momentum Tch as a function of the average to, for a few values 
of (T. At very low to, the coupled boost drops off rapidly because a 
large fraction of sightlines are optically thin - but in this case, the 
"boost" is negligible in either case (with or without leakage). At 
low to moderate to, the effective Tcff rises with r in a linear fash- 
ion as we would naively expect. At high-ro, however, the behavior 
depends on a, with a critical change in behavior around a = 1/2. 
This can be understood from the form of dP{T j to). 

For the distribution in Equation |B1[ the fraction of optically- 
thin sightlines (r < 1) scales simply as 

/.hm=^r-'/" (B2) 

So the average number of scatterings needed before a photon will 
"find" an optically thin sightline - and therefore have a high proba- 
bility of escaping the cloud - is, crudely, Nthm ~ 1 //thin oc rj^"^. On 
the other hand, the number of scatterings needed before the photon 
would diffuse out of the cloud assuming it did not find an optically 
thin sightline is just A'diffuse ~ tq. What matters for the behavior at 
high-ro is just this power-law falloff (exactly how we model the 
"core" and high-r part of the PDF make almost no difference to the 
values of Teff plotted in Figure [BT] at high-ro). 

If a < 1 /2, then, A'chin grows more rapidly than A'diffuse; in other 
words, the number of optically-thin sightlines falls off sufficiently 
rapidly at large ro that the typical photon will undergo its expected 
number of ^ To scatterings to diffuse out before it "leaks" out of an 
optically thin sightline, so that the coupled momentum T^ff oc ^/N oc 
rp. There is a linear normalization correction to rcff, which we can 
estimate analytically by considering the average distance travelled 
between scatterings, i.e. (A£) oc (r^') - this accounts for the fact 
that, on average, slightly more thin or thick sightlines allow for 
fewer or more scatterings before escape: we obtain the result that, 
for ro > 1, reff (1 - cr^) ro. 

If cr > 1/2, however, then A'diffuse grows more rapidly than 
Mhin; so the average photon will "leak out" after just Mhin scat- 
terings, before it can couple the full ro ~ VHofUix "boost." The 
actual coupled momentum should instead scale as rcff ~ x/AWi, 
which gives rcff oc t^^"^ . Again, we can analytically estimate the 
pre-factor for ro ^ 1, and obtain rcff — >■ a/2r[l/a] t^^"^ . The 
important point is that in this regime, the scaling is sub-linear in ro. 
There is still an approximately linear regime at moderate ro, but for 
very high ro, the "boost" becomes more limited. 

We stress that this behavior arises because the assumed 



dP{T \ To) behaves as a power-law at low r /ro - in other words, this 
allows for essentially the maximal "long tail" of low-r sightlines 
towards a central source with high average optical depth. Since 
iiP(dr/d^ I to/£q) reflects the local three-dimensional density dis- 
tribution, it might for example be more natural to assume it should 
have a lognormal form: 

dPiT\To)^^..p(-'^^)'^ (B3) 

The results of assuming this distribution are shown in Figure [B2] 
Given this distribution, the high-ro limit essentially always gives 
a linear scaling rcff oc ro. The reason is obvious given the ar- 
guments above - for a lognormal, the fraction of low-r sight- 
lines will decline much faster than a power-law, so independent 
of a, the probability of finding optically thin sightlines at r <C ro 
will fall much faster than rg^^. For a < 1, then, rcff — >■ ro when 
ro ^ 1. At very large a or more moderate ro, of course, r ~ 1 
may still fall within the "core" of the lognormal. For Equation |B 3 1 
the fraction of optically thin sightlines when ro ^ 1 is /thin — 
{a/V2^) In-' (ro) r-C"^"'/'^^'', or /,«„ oc r-C"^"'/'^'^''. The re- 
quirement that this fraction drop faster than r^-^ (to give rcff ~ ro) 
therfore becomes ro > exp (4cr"). For the reasonable values of a up 
to order unity, this is easily satisfied for high-ro. But if a were very 
large (say ~ 2), this rapidly becomes extremely large, and so we 
return to the Mwn < Miiffii.sc limit, and obtain the sub-linear scaling 
rcff ^ (27r)''''*[ln(ro)/(j]'/V(5'""°''''*"'. Even in this regime, how- 
ever, it is worth noting that at very high ro, the power-law model of 
Equation |Bl| still has a larger optically thin fraction - for a fixed a, 
that model gives a maximal effect of leakage. 

We can test the effects of this in our simulations by replac- 
ing the standard boost of ro with one of the appropriate rcff cal- 
culated above with some fixed a (using the curves in Fig |BI| to 
define an interpolation table). Obviously, for either distribution, a 
value of a < 0.5 will make no difference to any of our conclusions 
because rcff oc ro: the normalization correction is completely equiv- 
alent to variations in rip, discussed in the text (and in the regime of 
very low r, the boost scaling is not important). And for a lognor- 
mal distribution, any ct < I will yield identical results. We there- 
fore consider experiments with the exponential/power-law distri- 
bution and assumed a — 0.5, I.O and lognormal distribution with 
assumed a — 1.0, 2.0. For the power-law distribution with a — 0.5 
or lognormal with a — 1.0, rcff begins to deviate from ro at ro 3> 1, 
but the differences are sufficiently small that we do not see a large 
effect (they are roughly comparable to choices of 77,, — 2/3 and 
= 1 /2, respectively, and so change the expected SFRs only at the 
20 — 30% level). For the very large choices of a, however, we do 
expect and see some deviations. The equilibrium SFR is systemat- 
ically larger by a factor of ~ 2, similar to a small r^p ~ 1/4 (since 
the median riR ~ 30 — 50 becomes rcff ~ 10 here, this is expected). 
Visual inspection in these cases also confirms there are some small 
sub-regions in the galaxy nucleus where the gas consumption is 
near-runaway (these do not contain much of the mass, but they have 
the highest densities, riR ~ 100). Finally, we have also considered 
runs in which the power-law model is adopted, but with a taken 
from the local gas properties. Specifically, we take all of the gas 
inside the identified clump radius iJcinmp, and (knowing the density 
of each particle) compute the dispersion in In (p) which we use as 
a. We also add in quadrature a minimum a = 0.25 (0. 1 dex) which 
is about the minimum dispersion we see in ultra-high resolution 
simulations (in order to again be conservative and allow for signifi- 
cant leakage). The results of this run are quite similar to our default 
reft = ro model and/or the a = 0.5 model, which we expect since, as 
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noted in the text, the typical a we measure in simulation clumps is 
about 0.5. It is worth noting though, that the typical cr increases as 
we consider older and older stars, as a consequence of feedback in 
earlier stages driving out gas and "punching holes" in the gas distri- 
bution. Of course, the mean optical depth is also going down here, 
and we saw in the text that stars with the youngest ages < 1 Myr 
dominate the radiative momentum input. So accounting for leak- 
age accelerates the rate at which old stars luminosity can escape 
without coupling in the IR, but does not change our conclusions. 

We should also emphasize that other sub-grid effects could in 
fact raise Teff at fixed ro. This includes some non-trivial geometric 
cases where photons can be more efficiently trapped. Also, recall 
that we define to in the model as the globally averaged r out to 
a given radius oc Ms„cR^^, rather than the line-of-sight integrated 
r. If the gas within R is distributed with any average density profile 
that rises towards the stars, the appropriate Teff should be larger than 



that in Fig. BI (for a pure power-law profile p oc r^", this gives a 
factor oc 1/(1 — a) which is actually divergent for a > 1 ). In fact, if 
we calculate the true median line-of-sight integrated r in our ultra- 
high resolution simulations and compare it to the adopted to, the 
typical correction would amount to a "boost" of rjp ~ 2. And if we 
under-resolve collapse such that the gas "should" collapse a factor 
~ ip further in radius than our resolution limit allows, then Teff ~ 
tp^ To would be appropriate. It is difficult, therefore, to identify a 
"more accurate" model than our Teff oc to that is robust at the factor 
~ 2 level. 

We therefore conclude that photon leakage is unlikely to quali- 
tatively change our conclusions, given observationally and theoret- 
ically realistic distributions of column densities towards optically 
thick sources. However, it might be important in the most dense 
systems observed: starburst nuclei and AGN. The average IR op- 
tical depths in these regions can reach values > 100. The absolute 
value of the correction to Teff here could therefore be quite large - a 
factor ~ 10 rather than ~ 2, if cr is sufficiently large. Moreover, the 
sub-linear behavior of Teff could be very important, because these 
regions both have high-To and have dynamical times that are short 
relative to the stellar evolution timescale. In this joint limit, the lu- 
minosity required to support the system is TeffL/M oc (M/R^) or 
L/M oc To/Teff. When we are in the linear regime (Teff oc to), either 
because of low a or low to < 10, this implies that the system can 
self-regulate on both small and large scales once a fraction (a few 
percent) of the mass becomes stars. However, if Teff is significantly 
sub-linear, than the L/M needed to stabilize is a rising function of 
T - in other words, the system is vulnerable to runaway collapse 
( [Fall et al.|2010[ l. Such collapse could be quite interesting in these 
regions, however, since it would proceed with regions above a criti- 
cal To running away to turn entirely into stars, while neighboring re- 
gions that had smaller to do not collapse - and once the global L/M 
reached a given threshold, the low-density regions would be self- 
regulated. One might imagine a regime of global self-regulation 
on these scales, but without local self-regulation, perhaps making 
these regimes particularly interesting for the formation of globular 
clusters, dwarf galaxy nuclei, and/or super star clusters. 



